In [1]:
x = 3
y = 1.2
z = "Hi!"
Out[1]:
"Hi!"
In [2]:
typeof(x)
Out[2]:
Int64
In [3]:
typeof(y)
Out[3]:
Float64
In [4]:
typeof(z)
Out[4]:
String
In [5]:
x + 1
Out[5]:
4
In [6]:
# \mu + tab
μ = 0.0
σ = 2.0
xᵢ = 0
σ² = σ^2
Out[6]:
4.0
In [7]:
x = 1.0
y = 2.0
if x < y
    println("x is less than y.")
elseif x == y
    println("x is equal to y.")
else
    println("x is greter than y.")
end
x is less than y.
In [8]:
x = -3
# x < 0が成り立てば:の手前をyに代入する。成り立たねければ:の後をyに代入する
y = x < 0 ? "T" : "F"
println(y)
T
In [9]:
a = 0
b = 0
true && (a = 1) # (A && B)AがtrueならばBを実行
true || (b = 1) # (A || B)AがfalseならばBを実行
println("a = $(a), b = $(b)")
a = 1, b = 0
In [10]:
for i in 1:10
    println(i)
end
1
2
3
4
5
6
7
8
9
10
In [11]:
myinv(x) = 1 / x
Out[11]:
myinv (generic function with 1 method)
In [12]:
myinv(3)
Out[12]:
0.3333333333333333
In [13]:
function mymean(x, y)
    return (x + y) / 2
end
Out[13]:
mymean (generic function with 1 method)
In [14]:
mymean(1.0, 2.0)
Out[14]:
1.5
In [15]:
# 縦ベクトル
a = [1, 2, 3]
Out[15]:
3-element Vector{Int64}:
 1
 2
 3
In [16]:
# 横ベクトル
b = [1 2 3]
Out[16]:
1×3 Matrix{Int64}:
 1  2  3
In [17]:
c = Array{Float64}(undef, 3)
Out[17]:
3-element Vector{Float64}:
 2.2439480934e-314
 2.263943333e-314
 2.2426095865e-314
In [18]:
d = zeros(3)
Out[18]:
3-element Vector{Float64}:
 0.0
 0.0
 0.0
In [19]:
e = ones(3)
Out[19]:
3-element Vector{Float64}:
 1.0
 1.0
 1.0
In [20]:
# 0~1一様
f = rand(3)
Out[20]:
3-element Vector{Float64}:
 0.31385851194472925
 0.9294287946579434
 0.5799688477670801
In [21]:
# 標準正規分布
g = randn(3)
Out[21]:
3-element Vector{Float64}:
 -2.1397517175103573
 -0.3955230987292779
  0.31163264486238834
In [22]:
A = [1 2 3 4;
        5 6 7 8]
Out[22]:
2×4 Matrix{Int64}:
 1  2  3  4
 5  6  7  8
In [23]:
B = ones(2,4)
Out[23]:
2×4 Matrix{Float64}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
In [24]:
size(A)
Out[24]:
(2, 4)
In [25]:
length(A)
Out[25]:
8
In [26]:
A[2,1]
Out[26]:
5
In [27]:
A[:,1]
Out[27]:
2-element Vector{Int64}:
 1
 5
In [28]:
A[:,1:3]
Out[28]:
2×3 Matrix{Int64}:
 1  2  3
 5  6  7
In [29]:
[2*i for i in 1:5]
Out[29]:
5-element Vector{Int64}:
  2
  4
  6
  8
 10
In [30]:
[i*j for i in 1:3, j in 1:4]
Out[30]:
3×4 Matrix{Int64}:
 1  2  3   4
 2  4  6   8
 3  6  9  12
In [31]:
params1 = (1, 2, 3)
Out[31]:
(1, 2, 3)
In [32]:
f1(a, b, c) = a + b + c
f1(params1...)
Out[32]:
6
In [33]:
params2 = (1, 2, 3, 4)
Out[33]:
(1, 2, 3, 4)
In [34]:
f2(a, b, c, d) = a + b + c + d
f2(params2...)
Out[34]:
10
In [35]:
a = [1, 2, 3]
a + 1
MethodError: no method matching +(::Vector{Int64}, ::Int64)
For element-wise addition, use broadcasting with dot syntax: array .+ scalar
Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...) at operators.jl:560
  +(::T, ::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} at int.jl:87
  +(::Base.TwicePrecision, ::Number) at twiceprecision.jl:267
  ...

Stacktrace:
 [1] top-level scope
   @ In[35]:2
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [36]:
a .+ 1
Out[36]:
3-element Vector{Int64}:
 2
 3
 4
In [37]:
# ブロードキャスト
function add_one(x)
    x + 1
end
add_one(1.0)
Out[37]:
2.0
In [38]:
add_one.([1,2,3])
Out[38]:
3-element Vector{Int64}:
 2
 3
 4
In [39]:
x -> x + 1
Out[39]:
#5 (generic function with 1 method)
In [40]:
map(x -> x+1, [1,2,3])
Out[40]:
3-element Vector{Int64}:
 2
 3
 4
In [41]:
function test(maxiter)
    a = []
    for i in 1:maxiter
        push!(a, randn())
    end
    sum(a)
end
@time test(100000)
  0.008677 seconds (201.50 k allocations: 5.155 MiB, 36.59% compilation time)
Out[41]:
165.426717533415
In [42]:
mean([1,2,3,4,5])
UndefVarError: mean not defined

Stacktrace:
 [1] top-level scope
   @ In[42]:1
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [43]:
"""
よく使うパッケージ
Distirbutions.jl
ForwordDiff.jl
IJulia.jl
PyPlot.jl
"""
using Pkg
Pkg.add("Statistics")
cannot document the following expression:

using Pkg


Stacktrace:
 [1] error(::String, ::String)
   @ Base ./error.jl:42
 [2] top-level scope
   @ In[43]:1
 [3] eval
   @ ./boot.jl:360 [inlined]
 [4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [44]:
using Statistics
mean([1,2,3,4,5])
Out[44]:
3.0
In [45]:
Pkg.add("PyPlot")
UndefVarError: Pkg not defined

Stacktrace:
 [1] top-level scope
   @ In[45]:1
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [46]:
using PyPlot
In [47]:
# 1つ目のデータ
Y1 = [1, 7, 11, 13, 15, 16]

# 2つ目のデータ
Y2 = [15, 3, 13, 2, 7, 1]

# グラフを作成
fig, ax = subplots()

# 折れ線グラフを描画
ax.plot(Y1)
ax.plot(Y2)

# 軸の名称を設定
ax.set_xlabel("index")
ax.set_ylabel("Y")

# グラフのタイトルを設定
ax.set_title("my graph")
Out[47]:
PyObject Text(0.5, 1.0, 'my graph')
In [48]:
fig, axes = subplots(1, 2, figsize=(10,3))
axes[1].plot(Y1)
axes[2].plot(Y2)
Out[48]:
1-element Vector{PyCall.PyObject}:
 PyObject <matplotlib.lines.Line2D object at 0x7fe074f915e0>
In [49]:
# 2次関数の定義
quadraticf(x) = x^2

# -3から3まで、10個の等間隔の点列を生成
xs = range(-3, 3, length=10)

fig, ax = subplots()

# 関数の描画
ax.plot(xs, quadraticf.(xs), "-")

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("y=x^2")
Out[49]:
PyObject Text(0.5, 1.0, 'y=x^2')
In [50]:
# -3から3まで、10個の等間隔の点列を生成
xs = range(-3, 3, length=100)

fig, ax = subplots()

# 関数の描画
ax.plot(xs, quadraticf.(xs), "-")

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("y=x^2")
Out[50]:
PyObject Text(0.5, 1.0, 'y=x^2')
In [51]:
r = 1.0
fx(θ) = r*(θ-sin(θ))
fy(θ) = r*(1-cos(θ))
θs = range(-3pi, 3pi, length=100)
fig, ax = subplots()
ax.plot(fx.(θs), fy.(θs))
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("cycloid")
Out[51]:
PyObject Text(0.5, 1.0, 'cycloid')
In [52]:
# シグモイド関数を定義
sig(x) = 1 / (1 + exp(-x))

xs = range(-5, 5, length=100)
fig, ax = subplots()
ax.plot(xs, sig.(xs))
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("sigmoid")
Out[52]:
PyObject Text(0.5, 1.0, 'sigmoid')
In [53]:
# 可視化に使う乱数を生成
X = randn(1000)
Y = randn(1000)

fig, ax =subplots()

# ax.plot(X, Y, "o")などでも可
ax.scatter(X, Y)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("scatter")
Out[53]:
PyObject Text(0.5, 1.0, 'scatter')
In [54]:
X = randn(1000)
fig, ax = subplots()
ax.hist(X)
Out[54]:
([2.0, 22.0, 47.0, 145.0, 212.0, 246.0, 186.0, 86.0, 39.0, 15.0], [-3.3668178997339817, -2.7338823046905456, -2.1009467096471095, -1.468011114603673, -0.8350755195602368, -0.20213992451680074, 0.4307956705266358, 1.0637312655700715, 1.696666860613508, 2.3296024556569446, 2.96253805070038], (PyObject <matplotlib.patches.Rectangle object at 0x7fe05f453af0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f444b80>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f453fa0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f4641f0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f464400>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f464610>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f464820>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f464a30>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f464c40>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f464e50>))
In [55]:
X = randn(1000)
fig, ax = subplots()
ax.hist(X, bins=50)
Out[55]:
([3.0, 0.0, 2.0, 1.0, 0.0, 3.0, 4.0, 4.0, 5.0, 7.0  …  3.0, 5.0, 2.0, 1.0, 1.0, 3.0, 0.0, 0.0, 0.0, 1.0], [-3.3940216424986653, -3.2567722674011526, -3.1195228923036398, -2.9822735172061265, -2.8450241421086138, -2.707774767011101, -2.5705253919135878, -2.433276016816075, -2.296026641718562, -2.1587772666210494  …  2.233202736499363, 2.370452111596876, 2.5077014866943887, 2.6449508617919024, 2.782200236889415, 2.919449611986928, 3.0566989870844408, 3.1939483621819535, 3.3311977372794663, 3.46844711237698], (PyObject <matplotlib.patches.Rectangle object at 0x7fe05f60deb0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f4dff40>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61c3a0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61c5b0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61c7c0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61c9d0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61cbe0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61cdf0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61cfa0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62d250>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62d460>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62d670>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62d880>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62da90>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62dca0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62deb0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63b100>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63b310>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63b520>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63b730>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63b940>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63bb50>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63bd60>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63bf70>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64b1c0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64b3d0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64b5e0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64b7f0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64ba00>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64bc10>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64be20>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64bfd0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65b280>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65b490>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65b6a0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65b8b0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65bac0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65bcd0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65bee0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668130>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668340>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668550>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668760>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668970>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668b80>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668d90>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668fa0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f6791f0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f679400>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f679610>))
In [56]:
# 2変数関数を定義
fz(x, y) = exp(-(2x^2+y^2+x*y))

# x軸とy軸の可視化範囲を定義
xs = range(-1,1, length=100)
ys = range(-2,2, length=100)

fig, axes = subplots(1,2, figsize=(8,4))

# 1つ目のsubplotには数値をつける
cs = axes[1].contour(xs, ys, fz.(xs', ys))
axes[1].clabel(cs, inline=true)

# 2つ目のsubplotにはカラーバーをつける
cs = axes[2].contourf(xs, ys, fz.(xs', ys))
fig.colorbar(cs)

# subplots間の余白の大きさを自動調整
tight_layout()
In [57]:
a = [1, 2, 3]
2*a
Out[57]:
3-element Vector{Int64}:
 2
 4
 6
In [58]:
b = [4, 5, 6]
a + b
Out[58]:
3-element Vector{Int64}:
 5
 7
 9
In [59]:
a .* b
Out[59]:
3-element Vector{Int64}:
  4
 10
 18
In [60]:
#内積
sum(a .* b)
Out[60]:
32
In [61]:
a' * b
Out[61]:
32
In [62]:
A = [1 2 3;
        4 5 6]
B = [10 20 30;
        40 50 60]
A + B
Out[62]:
2×3 Matrix{Int64}:
 11  22  33
 44  55  66
In [63]:
A = [1 2;
        3 4;
        5 6]
B = [10 20 30 40;
        50 60 70 80]
C = A * B
Out[63]:
3×4 Matrix{Int64}:
 110  140  170  200
 230  300  370  440
 350  460  570  680
In [64]:
M = size(A,1) #Aの行数
N = size(B, 2) #Bの列数

# M * Nの行列を作成
C = [sum(A[i,:] .* B[:,j]) for i in 1:M, j in 1:N]
Out[64]:
3×4 Matrix{Int64}:
 110  140  170  200
 230  300  370  440
 350  460  570  680
In [65]:
B * A
DimensionMismatch("matrix A has dimensions (2,4), matrix B has dimensions (3,2)")

Stacktrace:
 [1] _generic_matmatmul!(C::Matrix{Int64}, tA::Char, tB::Char, A::Matrix{Int64}, B::Matrix{Int64}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
   @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:814
 [2] generic_matmatmul!(C::Matrix{Int64}, tA::Char, tB::Char, A::Matrix{Int64}, B::Matrix{Int64}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
   @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:802
 [3] mul!
   @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:302 [inlined]
 [4] mul!
   @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:275 [inlined]
 [5] *(A::Matrix{Int64}, B::Matrix{Int64})
   @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:153
 [6] top-level scope
   @ In[65]:1
 [7] eval
   @ ./boot.jl:360 [inlined]
 [8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [66]:
A = [1 2;
        3 4;
        5 6]
I = [1 0 0;
        0 1 0;
        0 0 1]
I * A
Out[66]:
3×2 Matrix{Int64}:
 1  2
 3  4
 5  6
In [67]:
A = [1 2 3;
        4 5 6]
A'
Out[67]:
3×2 adjoint(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5
 3  6
In [68]:
a = [1, 2, 3]
b = [5, 7]
a * b'
Out[68]:
3×2 Matrix{Int64}:
  5   7
 10  14
 15  21
In [69]:
f2(x, y) = 2*x + y
f2.(a, b')
Out[69]:
3×2 Matrix{Int64}:
  7   9
  9  11
 11  13
In [70]:
A = [1 2;
        3 4]
B = inv(A)
Out[70]:
2×2 Matrix{Float64}:
 -2.0   1.0
  1.5  -0.5
In [71]:
A * B
Out[71]:
2×2 Matrix{Float64}:
 1.0          0.0
 8.88178e-16  1.0
In [72]:
B * A
Out[72]:
2×2 Matrix{Float64}:
 1.0          0.0
 2.22045e-16  1.0
In [73]:
println(B)
[-1.9999999999999996 0.9999999999999998; 1.4999999999999998 -0.4999999999999999]
In [74]:
A = Rational{Int}[1 2;
                                3 4]
B = inv(A)
println(A * B)
println(B * A)
Rational{Int64}[1//1 0//1; 0//1 1//1]
Rational{Int64}[1//1 0//1; 0//1 1//1]
In [75]:
A = Rational{Int}[1 2;
                                3 4]
sol = inv(A) * [-1, 1]
Out[75]:
2-element Vector{Rational{Int64}}:
  3//1
 -2//1
In [76]:
using Statistics
In [77]:
X = rand(5)
Out[77]:
5-element Vector{Float64}:
 0.19604681062925633
 0.15039929085606318
 0.6494316784702181
 0.5041200110315147
 0.2762255348491489
In [78]:
Y = rand(2, 5)
Out[78]:
2×5 Matrix{Float64}:
 0.880842  0.610478  0.409949  0.1643    0.516474
 0.488594  0.386113  0.834355  0.291409  0.122172
In [79]:
println(sum(X))
println(mean(X))
1.7762233258362012
0.35524466516724024
In [80]:
println(sum(Y))
println(sum(Y, dims=1))
println(sum(Y, dims=2))
4.704685348250634
[1.3694354496341408 0.996590834143914 1.2443035609926136 0.4557090946317597 0.6386464088482056]
[2.582042658726274; 2.1226426895243593]
In [81]:
println(std(X))
println(std(X)^2)
println(var(X))
0.21345929039558376
0.045564868656186155
0.045564868656186155
In [82]:
cov(Y,dims=1)
Out[82]:
5×5 Matrix{Float64}:
  0.0769292   0.0440033  -0.0832362  -0.024929     0.0773321
  0.0440033   0.0251697  -0.0476108  -0.0142593    0.0442337
 -0.0832362  -0.0476108   0.0900602   0.0269728   -0.0836721
 -0.024929   -0.0142593   0.0269728   0.00807828  -0.0250596
  0.0773321   0.0442337  -0.0836721  -0.0250596    0.077737
In [83]:
cov(Y,dims=2)
Out[83]:
2×2 Matrix{Float64}:
 0.0692436   0.00573912
 0.00573912  0.0706695
In [84]:
Pkg.add("Distributions")
UndefVarError: Pkg not defined

Stacktrace:
 [1] top-level scope
   @ In[84]:1
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [85]:
using Distributions
μ = 1.5
σ = 2.0
Z = rand(Normal(μ, σ), 10000)
Z
Out[85]:
10000-element Vector{Float64}:
 -2.878438829554005
  1.9253341325276128
  0.7912537284479731
 -0.28884382408978526
  7.4798218320469
  1.4923605551865857
  3.0248292221970154
  1.039516574189463
  5.394882281609656
  0.14356031974931138
  0.4469234226249874
  1.8896840550067933
  3.5436087455003693
  ⋮
 -0.14534552137994927
 -1.7913541860662607
  2.4454416739331584
  0.007851113359583328
 -0.6249474431000608
  2.12050537813308
 -1.0890719963791056
 -1.279649551121989
  0.9741969853930287
  2.583001380041292
  1.5100434409297825
  0.34536356361712217
In [86]:
println(mean(Z))
println(std(Z))
1.524948267501824
1.9954566091407018
In [87]:
Pkg.add("PyPlot")
UndefVarError: Pkg not defined

Stacktrace:
 [1] top-level scope
   @ In[87]:1
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [88]:
using PyPlot

# quadraticf(x)を2次関数として定義
quadraticf(x) = -(x + 1)*(x-1)

# hを微小な量として定義(10のマイナス10乗)
h = 1.0e-10

# 導関数quadraticf'の近似式
quadraticff(a) = (quadraticf(a+h) - quadraticf(a)) / h

# 関数の可視化範囲
xs = range(-1, 1, length=100)

fig, axes = subplots(2,1, figsize=(4,6))

# 関数のプロット
axes[1].plot(xs, quadraticf.(xs), "b")
axes[1].grid()
axes[1].set_xlabel("x"), axes[1].set_ylabel("y")
axes[1].set_title("function f(x)")

# 導関数のプロット
axes[2].plot(xs, quadraticff.(xs), "r")
axes[2].grid()
axes[2].set_xlabel("x"), axes[2].set_ylabel("y")
axes[2].set_title("derivative f'(x)")

tight_layout()
In [89]:
# グラフを可視化する解像度
L = 10

# f2(x)を可視化する範囲
xs1 = range(-1, 1, length=L)
xs2 = range(-1, 1, length=L)

# 2変数関数の定義
f2(x) = -(x .+ 1)' * (x .- 1)

# 2変数関数の勾配
∇f2(x) = -2x

fig, axes = subplots(1,2, figsize=(8,4))

# 関数の等高線の可視化
cs = axes[1].contour(xs1, xs2, [f2([x1,x2]) for x1 in xs1, x2 in xs2]')
axes[1].clabel(cs, inline=true)
axes[1].set_xlabel("x1"), axes[1].set_ylabel("x2")
axes[1].set_title("f2(x)")

# 勾配ベクトルの計算と可視化
vec1 = [∇f2([x1,x2])[1] for x1 in xs1, x2 in xs2]
vec2 = [∇f2([x1,x2])[2] for x1 in xs1, x2 in xs2]
axes[2].quiver(repeat(xs1, 1, L), repeat(xs2, L, 1), vec1, vec2)
axes[2].set_xlabel("x1"), axes[2].set_ylabel("x2")
axes[2].set_title("∇f2(x)")

tight_layout()
In [90]:
Pkg.add("ForwardDiff")
UndefVarError: Pkg not defined

Stacktrace:
 [1] top-level scope
   @ In[90]:1
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [91]:
using ForwardDiff
In [92]:
# 2変数関数の定義
f3(x) = -(x + 1) * (x - 1)

# 自動微分によって導関数f'(x)を求める
f33(x) = ForwardDiff.derivative(f3, x)

# 関数の可視化範囲
xs = range(-1, 1, length=100)

fig, axes = subplots(2,1, figsize=(4,6))

# 関数のプロット
axes[1].plot(xs, f3.(xs), "b")
axes[1].grid()
axes[1].set_xlabel("x"), axes[1].set_ylabel("y")
axes[1].set_title("functin f(x)")

# 導関数のプロット
axes[2].plot(xs, f33.(xs), "r")
axes[2].grid()
axes[2].set_xlabel("x"), axes[1].set_ylabel("y")
axes[2].set_title("derivative f'(x)")

tight_layout()
In [93]:
fig, ax = subplots()
xs = range(0, 2pi*3, length=100)

# sin(x)をプロット
ax.plot(xs, sin.(xs), color="b", label="sin(x)")

# 導関数をプロット
ax.plot(xs, map(x -> ForwardDiff.derivative(sin, x), xs), color="r", label="sin'(x)")
ax.grid()
ax.set_xlabel("x"), ax.set_ylabel("y")
ax.legend()
Out[93]:
PyObject <matplotlib.legend.Legend object at 0x7fe0608ff310>
In [94]:
# シグモイド関数を定義
sig(x) = 1/(1 + exp(-x))

xs = range(-5, 5, length=100)
fig, axes = subplots(2, 1, figsize=(4,6))

# シグモイド関数をプロット
axes[1].plot(xs, sig.(xs), color="b")
axes[1].set_xlabel("x"), axes[1].set_ylabel("y")
axes[1].set_title("sig(x)")
axes[1].grid()

# 導関数をプロット
axes[2].plot(xs, map(x -> ForwardDiff.derivative(sig, x), xs), color="r")
axes[2].set_xlabel("x"), axes[2].set_ylabel("y")
axes[2].set_title("sig'(x)")
axes[2].grid()
In [95]:
# 最大値を探したい関数
x_opt = 0.50
quadraticf(x) = -2(x - x_opt)^2

fig, ax = subplots()
xs = range(-3, 3, length=100)
ax.plot(xs, quadraticf.(xs), label="function")
ax.plot(x_opt, quadraticf(x_opt), "rx", label="optimal")
ax.set_xlabel("x"), ax.set_ylabel("y")
ax.grid()
ax.legend()
Out[95]:
PyObject <matplotlib.legend.Legend object at 0x7fe060e2c3a0>
In [96]:
# 1変数関数の最適化
function gradient_method_1dim(f, x_init, η, maxiter)
    # 最適化過程のパラメータを格納する配列
    x_seq = Array{typeof(x_init), 1}(undef, maxiter)
    
    # 勾配
    ff(x) = ForwardDiff.derivative(f, x)
    
    # 初期値
    x_seq[1] = x_init
    
    # メインの最適解ループ
    for i in 2:maxiter
        x_seq[i] = x_seq[i-1] + η*ff(x_seq[i-1])
    end
    
    x_seq

end
Out[96]:
gradient_method_1dim (generic function with 1 method)
In [97]:
# 探索範囲の初期値
x_init = -2.5

# 探索の繰り返し数
maxiter = 20

# ステップサイズ
η = 0.1

# 最適化計算を実行
x_seq = gradient_method_1dim(quadraticf, x_init, η, maxiter)
quadraticf_seq = quadraticf.(x_seq)
Out[97]:
20-element Vector{Float64}:
 -18.0
  -6.479999999999999
  -2.3327999999999993
  -0.8398079999999998
  -0.30233087999999986
  -0.10883911679999994
  -0.03918208204799999
  -0.014105549537279988
  -0.0050779978334207915
  -0.0018280792200314835
  -0.0006581085192113332
  -0.0002369190669160795
  -8.529086408978862e-5
  -3.07047110723239e-5
  -1.1053695986036396e-5
  -3.979330554973165e-6
  -1.4325589997903394e-6
  -5.157212399245448e-7
  -1.8565964637284966e-7
  -6.683747269421776e-8
In [98]:
# 目的関数の値をステップごとにプロット
fig, ax = subplots(figsize=(8,3))
ax.plot(quadraticf_seq)
ax.set_xlabel("iteration"), ax.set_ylabel("f")
ax.grid()

fig,axes = subplots(1,2,figsize=(8,3))

# 関数のプロット
axes[1].plot(xs, quadraticf.(xs), label="function")

# 探索の過程
axes[1].plot(x_seq, quadraticf.(x_seq), "bx", label="sequence")

# 最適値
axes[1].plot(x_opt, quadraticf(x_opt), "rx", label="optimal")
axes[1].set_xlabel("x"), axes[1].set_ylabel("y")
axes[1].grid()
axes[1].legend()

# 探索の過程をステップごとにプロット
axes[2].plot(1:maxiter, x_seq, "bx-", label="sequence")
axes[2].hlines(x_opt, 0, maxiter, ls="--", label="x_opt")
axes[2].set_xlim([0, maxiter])
axes[2].set_xlabel("iteration"), axes[2].set_ylabel("x")
axes[2].grid()
axes[2].legend()

tight_layout()
In [99]:
# 2次関数を定義
x_opt = [0.50, 0.25]
f2(x) = -sqrt(0.05 + (x[1] - x_opt[1])^2) - (x[2] - x_opt[2])^2

# 関数を等高線として可視化
L = 100
xs1 = range(-1, 1, length=L)
xs2 = range(-1, 1, length=L)
fig, ax = subplots()
ax.contour(xs1, xs2, [f2(x1, x2) for x1 in xs1, x2 in xs2]')
ax.plot(x_opt[1], x_opt[2], color="r", marker="x", label="optimal")
ax.set_xlabel("x1"), ax.set_ylabel("x2")
ax.grid()
ax.legend()
Out[99]:
PyObject <matplotlib.legend.Legend object at 0x7fe060f2dd00>
In [100]:
# 多変数関数のための勾配法
function gradient_method(f, x_init, η, maxiter)
    # 最適化過程のパラメータを格納する配列
    x_seq = Array{typeof(x_init[1]), 2}(undef, length(x_init), maxiter)
    
    # 勾配
    ∇f(x) = ForwardDiff.gradient(f, x)
    
    # 初期値
    x_seq[:, 1] = x_init
    
    # メインの最適化ループ
    for i in 2:maxiter
        x_seq[:, i] = x_seq[:, i-1] + η*∇f(x_seq[:, i-1])
    end
    
    x_seq
end
    
# パラメータの設定
x_init = [-0.75, -0,75]

maxiter = 20
η=0.1

# 最適化の実行
x_seq = gradient_method(f2, x_init, η, maxiter)
f_seq = [f2(x_seq[:,i]) for i in 1:maxiter]
Out[100]:
20-element Vector{Float64}:
 -1.3323425099200294
 -1.2130713587869222
 -1.10246741292214
 -0.9977733058718217
 -0.8973440357566991
 -0.8003014554196972
 -0.7063454734180662
 -0.6156878830524037
 -0.5291020764458453
 -0.448101220575972
 -0.3752360053913261
 -0.3143005037996588
 -0.26957995148648584
 -0.24287708480788592
 -0.23050261870859573
 -0.2258762089573429
 -0.22433899880772853
 -0.22384834870321288
 -0.22369132580962012
 -0.22363942129961714
In [101]:
# 目的関数の値をステップごとにプロット
fig, ax = subplots(figsize=(8,3))
ax.plot(f_seq)
ax.set_xlabel("iteration"), ax.set_ylabel("f")
ax.grid()

fig, axes = subplots(1, 2, figsize=(8, 3))

# 等高線図で関数を可視化
axes[1].contour(xs1, xs2, [f2(x1, x2) for x1 in xs1, x2 in xs2]')

# 最適化の過程
axes[1].plot(x_seq[1, :], x_seq[2,:], ls="--", marker="x", label="sequences")
axes[1].plot(x_opt[1], x_opt[2], color="r", marker="x", label="optimal")
axes[1].set_xlabel("x1"), axes[1].set_ylabel("x2")
axes[1].grid()
axes[1].legend()

# ステップごとの最適化の過程
axes[2].plot(1:maxiter, x_seq[1,:], color="b", marker="o", label="x[1]")
axes[2].plot(1:maxiter, x_seq[2,:], color="r", marker="^", label="x[2]")
axes[2].hlines(x_opt[1], 0, maxiter, color="b", ls="--", label="x_opt[1]")
axes[2].hlines(x_opt[2], 0, maxiter, color="r", ls="-", label="x_opt[2]")
axes[2].set_xlabel("iteration")
axes[2].set_ylabel("x1, x2")
axes[2].grid()
axes[2].legend()

tight_layout()
In [102]:
# 目的関数の定義
f_complex(x) = 0.3*cos(3pi*x) - x^2

# 最適化のパラメータ設定
maxiter = 20
η = 0.01

# 初期値a
x_init_a = -0.8
x_seq_a = gradient_method_1dim(f_complex, x_init_a, η, maxiter)
f_seq_a = f_complex.(x_seq_a)

# 初期値b
x_init_b = 0.25
x_seq_b = gradient_method_1dim(f_complex, x_init_b, η, maxiter)
f_seq_b = f_complex.(x_seq_b)

# 目的関数の値をステップごとにプロット
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(f_seq_a)
axes[1].set_xlabel("iteration"), axes[1].set_ylabel("f")
axes[1].set_title("f(x) (a)")
axes[1].grid()
axes[2].plot(f_seq_b)
axes[2].set_xlabel("iteration"), axes[2].set_ylabel("f")
axes[2].set_title("f(x) (b)")
axes[2].grid()
tight_layout()

# 関数を可視化する範囲
xs = range(-1, 1, length=100)

# 最適化の過程
fig, axes = subplots(1, 2, figsize=(12,4))
axes[1].plot(xs, f_complex.(xs), label="f_complex")
axes[1].plot(x_seq_a, f_complex.(x_seq_a), color="b", marker="x", label="x sequences")
axes[1].set_xlabel("x"), axes[1].set_ylabel("y")
axes[1].grid()
axes[1].set_title("initial value a")
axes[1].legend()
axes[2].plot(xs, f_complex.(xs), label="f_complex")
axes[2].plot(x_seq_b, f_complex.(x_seq_b), color="b", marker="x", label="x sequences")
axes[2].set_xlabel("x"), axes[2].set_ylabel("y")
axes[2].grid()
axes[2].set_title("initial value b")
axes[2].legend()

tight_layout()
In [103]:
# 学習用の入力値集合
X_obs = [1, 2, 4]

# 学習用の出力値集合
Y_obs = [10, 35, 72]

# データの可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
ax.set_xlabel("x"), ax.set_ylabel("y")
ax.set_title("data")
ax.grid()
In [104]:
# 適当な傾きパラメータおよび切片パラメータを設定
w = [12.0, 10.0]

# 予測に使う関数
ff(x) = w[1]*x + w[2]

# 関数を可視化する範囲
xs = range(0, 5, length=100)

# データと予測関数の可視化
fig, ax = subplots()
ax.plot(xs, ff.(xs), label="function")
ax.scatter(X_obs, Y_obs)
ax.set_xlabel("x"), ax.set_ylabel("y")
ax.grid()
ax.legend()
Out[104]:
PyObject <matplotlib.legend.Legend object at 0x7fe062cb7a00>
In [105]:
# 誤差関数の定義
E(w) = sum([(Y_obs[n] - (w[1]*X_obs[n] + w[2]))^2 for n in 1:length(X_obs)])
Out[105]:
E (generic function with 1 method)
In [106]:
# 最適化するパラメータの初期値
w_init = [0.0, 0.0]

# 最適化計算の回数
maxiter = 500

#学習率
η = 0.01

# 最適化の実行、最大化アルゴリズムなので、-Eを目的関数にする
F(w) = -E(w)
w_seq = gradient_method(F, w_init, η, maxiter)
f_seq = [F(w_seq[:,i]) for i in 1:maxiter]
Out[106]:
500-element Vector{Float64}:
 -6509.0
 -1939.5819999999997
  -644.3813987199998
  -275.7910161687038
  -169.46811483893487
  -137.41138630873562
  -126.41782512660872
  -121.4305759666068
  -118.18766202638302
  -115.48193910495918
  -112.97076516030776
  -110.556165388921
  -108.20942661889302
     ⋮
   -12.072442904302092
   -12.072418776699926
   -12.072395223013139
   -12.072372229590176
   -12.072349783104178
   -12.072327870545308
   -12.072306479213296
   -12.072285596709893
   -12.072265210931684
   -12.072245310063245
   -12.072225882570269
   -12.072206917192597
In [107]:
w1, w2 = w_seq[:, end]
println("w_1 = $(w1), w_2 = $(w2)")
w_1 = 20.345436819385633, w_2 = -8.465882327774901
In [108]:
# 目的関数の値をステップごとにプロット
fig, ax = subplots(figsize=(8,4))
ax.plot(f_seq)
ax.set_xlabel("iteration"), ax.set_ylabel("f")
ax.grid()

# 予測に使う関数
ff(x) = w1*x + w2

# 関数を可視化する範囲
xs = range(0, 5, length=100)

# データと予測関数をプロット
fig, ax = subplots()
ax.plot(xs, ff.(xs), label="f(x)")
ax.scatter(X_obs, Y_obs, label="data")
ax.set_xlabel("x"), ax.set_ylabel("y")
ax.grid()
ax.legend()
Out[108]:
PyObject <matplotlib.legend.Legend object at 0x7fe0615cef10>
In [109]:
function linear_fit(Y, X)
    N = length(Y)
    w1 = sum((Y .- mean(Y)) .* X) / sum((X .- mean(X)).*X)
    w2 = mean(Y) - w1*mean(X)
    w1, w2
end
w1, w2 = linear_fit(Y_obs, X_obs)
println("w1 = $(w1), w2 = $(w2)")
w1 = 20.35714285714286, w2 = -8.500000000000014
In [110]:
function approx_integration(x_range, f)
    # 幅
    Δ = x_range[2] - x_range[1]
    
    # 近似された面積と幅を返す
    sum([f(x) * Δ for x in x_range]), Δ
end
Out[110]:
approx_integration (generic function with 1 method)
In [111]:
# 面積を近似したい関数
ff(x) = exp(-x^2)

# 等間隔の配列を用意
#x_range = range(-2, 2, length=10)
x_range = range(-5, 5, length=100)

# 積分近似の実行
approx, Δ = approx_integration(x_range, ff)

# 近似値(apporox)と厳密値(exact)の比較
println("approx = $(approx)")
println("exact = $(sqrt(pi))")

# 関数と近似結果の可視化
fig, ax = subplots(figsize=(8,4))
xs = range(-5, 5, length=1000)
ax.plot(xs, ff.(xs), "r--")
for x in x_range
    ax.fill_between([x - Δ/2, x + Δ/2], [ff(x), ff(x)], [0, 0], facecolor="b", edgecolor="k", alpha=0.5)
end
# ax.set_xlabel("x"), ax.set_ylabel("y")
# ax.grid()
approx = 1.7724538509039651
exact = 1.7724538509055159
In [112]:
# 積分対象の2変数関数
D = 2
A = [0.5 0.3 
        0.3 1.0]
f2(x) = exp(-0.5*x'A*x)

# 20x20の区画に分割
L = 20
xs1 = range(-5, 5, length=L)
xs2 = range(-5, 5, length=L)

fig, axes = subplots(1, 2, figsize=(8,4))

# 等高線図で可視化
cs = axes[1].contour(repeat(xs1, 1, L), repeat(xs2', L, 1), [f2([x1, x2]) for x1 in xs1, x2 in xs2]')
axes[1].clabel(cs)
axes[1].grid()
axes[1].set_xlabel("x1"), axes[1].set_ylabel("x2")

# カラー表示
cs = axes[2].imshow([f2([x1, x2]) for x1 in xs1, x2 in xs2], origin="lower")
fig.colorbar(cs)
axes[2].set_xlabel("x1"), axes[2].set_ylabel("x2")

tight_layout()
In [113]:
function approx_integration_2dim(x_range, f)
    Δ = x_range[2] - x_range[1]
    sum([f([x1, x2]) * Δ^2 for x1 in x_range, x2 in x_range]), Δ
end
Out[113]:
approx_integration_2dim (generic function with 1 method)
In [114]:
# det関数を使うためにLinearAlgebraパッケージを使用
Pkg.add("LinearAlgebra")
using LinearAlgebra

# 20x20の区間に分割
L = 20
x_range = range(-5, 5, length=L)
approx, Δ = approx_integration_2dim(x_range, f2)

println("approx = $(approx)")
println("exact = $(sqrt((2*pi)^D/det(A)))")
UndefVarError: Pkg not defined

Stacktrace:
 [1] top-level scope
   @ In[114]:2
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [115]:
# 1000x1000の区間に分割
L = 1000
x_range = range(-100, 100, length=L)
approx, Δ = approx_integration_2dim(x_range, f2)

println("approx = $(approx)")
println("exact = $(sqrt((2*pi)^D/det(A)))")
approx = 9.812686860654521
UndefVarError: det not defined

Stacktrace:
 [1] top-level scope
   @ In[115]:7
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [116]:
using Distributions

# パラメータが0.5のベルヌーイ分布を定義
bern = Bernoulli(0.5)

# 乱数を10個生成
X = rand(bern, 10)
Out[116]:
10-element Vector{Bool}:
 1
 1
 1
 0
 1
 1
 0
 0
 1
 1
In [117]:
#パラメータを変更
bern = Bernoulli(0.9)

X = rand(bern, 10)
Out[117]:
10-element Vector{Bool}:
 1
 1
 1
 1
 1
 1
 0
 0
 0
 1
In [118]:
bag(x::Bool) = x == 1 ? "A" : "B"
ball(y::Bool) = y == 1 ? "赤" : "白"
X = bag.(rand(bern, 10))
Out[118]:
10-element Vector{String}:
 "A"
 "A"
 "A"
 "A"
 "A"
 "B"
 "A"
 "A"
 "A"
 "A"
In [119]:
function sample()
    # 袋 の選択はそれぞれ1/2の確率
    x = bag.(rand(Bernoulli(1//2)))
    
    # 袋がAであれば、赤玉が出る確率は1/5
    # 袋がBであれば、赤玉が出る確率は3/5
    μ = x=="A" ? 1//5 : 3//5
    
    # 玉の抽出
    y = ball.(rand(Bernoulli(μ)))
    
    x, μ, y
end
Out[119]:
sample (generic function with 1 method)
In [120]:
for _ in 1:10
    x, μ, y = sample()
    println("袋: $(x), 玉: $(y)")
end
袋: A, 玉: 白
袋: B, 玉: 赤
袋: A, 玉: 赤
袋: A, 玉: 赤
袋: B, 玉: 赤
袋: A, 玉: 赤
袋: B, 玉: 白
袋: B, 玉: 赤
袋: A, 玉: 白
袋: B, 玉: 赤
In [121]:
maxiter = 100
result = []
for _ in 1:maxiter
    x, μ, y = sample()
    push!(result, y)
end
mean(result .== "赤")
Out[121]:
0.45
In [122]:
maxiter = 1000000
result = []
for _ in 1:maxiter
    x, μ, y = sample()
    push!(result, y)
end
mean(result .== "赤")
Out[122]:
0.399211
In [123]:
# 観測値(赤玉)
y_obs = "赤"

maxiter = 1_000_000
x_results = []
for _ in 1:maxiter
    x, μ, y = sample()
    
    # 生成されたyが観測と一致する場合のみ追加
    y == y_obs && push!(x_results, x)
end

# 受容率(観測と一致した割合)
println("acceptance rate = $(length(x_results)/maxiter)")

# 玉が赤だった場合の袋の条件付き分布
println("p(x=A|y=赤) : approx = $(mean(x_results .== "A")) (exact=$(1/4))")
println("p(x=B|y=赤) : approx = $(mean(x_results .== "B")) (exact=$(3/4))")
acceptance rate = 0.399853
p(x=A|y=赤) : approx = 0.24998436925570147 (exact=0.25)
p(x=B|y=赤) : approx = 0.7500156307442986 (exact=0.75)
In [124]:
function sample(N)
    x = bag(rand(Bernoulli(1//2)))
    μ = x=="A" ? 1//5 : 3//5
    
    # 玉をN回復元抽出する
    Y = ball.(rand(Bernoulli(μ), N))
    
    x, μ, Y
end
Out[124]:
sample (generic function with 2 methods)
In [125]:
maxiter = 10_000
Y_obs = ["赤", "赤", "白"]
x_results = []
for _ in 1:maxiter
    x, μ, Y = sample(3)
    
    # 3つの玉が完全に一致する場合のみ受容
    Y == Y_obs && push!(x_results, x)
end
println("acceptance rate = $(length(x_results)/maxiter)")
println("p(x=A|y_1=赤, y_2=赤, y_3=白) : " *
            "approx = $(mean(x_results .== "A")) (exact=$(2/11))")
println("p(x=B|y_1=赤, y_2=赤, y_3=白) : " *
            "approx = $(mean(x_results .== "B")) (exact=$(9/11))")
acceptance rate = 0.09
p(x=A|y_1=赤, y_2=赤, y_3=白) : approx = 0.19444444444444445 (exact=0.18181818181818182)
p(x=B|y_1=赤, y_2=赤, y_3=白) : approx = 0.8055555555555556 (exact=0.8181818181818182)
In [126]:
maxiter = 10_000
x_results = []
for _ in 1:maxiter
    x, μ, Y = sample(3)
    
    # 赤玉の個数さえ一致すれば受容するように修正
    sum(Y.=="赤") == sum(Y_obs.=="赤") && push!(x_results, x)
end
println("acceptance rate = $(length(x_results)/maxiter)")
println("p(x=A|y_1=赤, y_2=赤, y_3=白) : " *
            "approx = $(mean(x_results .== "A")) (exact=$(2/11))")
println("p(x=B|y_1=赤, y_2=赤, y_3=白) : " *
            "approx = $(mean(x_results .== "B")) (exact=$(9/11))")
acceptance rate = 0.2638
p(x=A|y_1=赤, y_2=赤, y_3=白) : approx = 0.17930250189537528 (exact=0.18181818181818182)
p(x=B|y_1=赤, y_2=赤, y_3=白) : approx = 0.8206974981046247 (exact=0.8181818181818182)
In [127]:
using Distributions
using PyPlot
function set_options(ax, xlabel, ylabel, title;
                                        grid=true, gridy=false, legend=false)
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)
    ax.set_title(title)
    if grid
        if gridy
            ax.grid(axis="y")
        else
            ax.grid()
        end
    end
    legend && ax.legend()
end
Out[127]:
set_options (generic function with 1 method)
In [128]:
# 分布を作る
d = Bernoulli(0.3)
Out[128]:
Bernoulli{Float64}(p=0.3)
In [129]:
x = rand(d)
x
Out[129]:
true
In [130]:
println(x)
println(Int(x))
true
1
In [131]:
# 10個の独立なサンプルを得る
X = rand(d, 10)
X'
Out[131]:
1×10 adjoint(::Vector{Bool}) with eltype Bool:
 0  1  1  0  0  0  0  0  1  0
In [132]:
# 1が生成される確率
println(pdf(d, 1))

# 0が生成する確率
println(pdf(d, 0))

# -1が生成される確率(起こり得ないのでゼロになる)
println(pdf(d, -1))
0.3
0.7
0.0
In [133]:
fig, ax = subplots()
ax.bar([0, 1], pdf.(d, [0,1]))
set_options(ax, "x", "probability", "probability mass function";
                        gridy=true)
Out[133]:
false
In [134]:
println("mean = $(mean(d)), std = $(std(d))")
mean = 0.3, std = 0.458257569495584
In [135]:
μ = 0.3
println("std = $(sqrt(μ*(1-μ)))")
std = 0.458257569495584
In [136]:
X = rand(d, 10_000)
println("mean ∼ $(mean(X)), std ∼ $(std(X))")
mean ∼ 0.3019, std ∼ 0.4591050726650431
In [137]:
d = Binomial(20, 0.3)
Out[137]:
Binomial{Float64}(n=20, p=0.3)
In [138]:
x = rand(d)
Out[138]:
6
In [139]:
X = rand(d, 100)
fig, ax = subplots()
ax.hist(X)
set_options(ax, "x", "freqency", "histogram";
                        gridy=true)
Out[139]:
false
In [140]:
println("mean (exact) = $(mean(d)), mean (approx) = $(mean(X))")
mean (exact) = 6.0, mean (approx) = 6.04
In [141]:
println("std (exact) = $(std(d)), std (approx) = $(std(X))")
std (exact) = 2.0493901531919194, std (approx) = 1.9741767239469057
In [142]:
xs = range(0, 20, length=21)
fig, ax = subplots()
ax.bar(xs, pdf.(d, xs))
set_options(ax, "x", "freqency", "probability mass function";
                        gridy=true)
Out[142]:
false
In [143]:
# 平均パラメータのリスト
μs = [0.2, 0.4, 0.6, 0.8]

fig, axes = subplots(2, 2, sharey=true, figsize=(8,6))
for (i, ax) in enumerate(axes)
    μ = μs[i]
    d = Binomial(20, μ)
    ax.bar(xs, pdf.(d, xs))
    set_options(ax, "x", "probability", "μ=$(μ)";
                            gridy=true)
end
tight_layout()
In [144]:
M = 10
d = Multinomial(M, [0.5, 0.3, 0.2])
Out[144]:
Multinomial{Float64, Vector{Float64}}(n=10, p=[0.5, 0.3, 0.2])
In [145]:
x = rand(d)
x
Out[145]:
3-element Vector{Int64}:
 6
 0
 4
In [146]:
X = rand(d, 100)
X
Out[146]:
3×100 Matrix{Int64}:
 5  5  6  3  5  3  6  4  4  4  3  4  4  …  5  3  5  5  5  3  4  4  5  6  6  4
 5  3  4  4  4  2  3  5  2  2  3  3  5     5  4  3  4  3  6  3  2  3  3  3  1
 0  2  0  3  1  5  1  1  4  4  4  3  1     0  3  2  1  2  1  3  4  2  1  1  5
In [147]:
mean(X, dims=2)
Out[147]:
3×1 Matrix{Float64}:
 4.59
 3.45
 1.96
In [148]:
mean(d)
Out[148]:
3-element Vector{Float64}:
 5.0
 3.0
 2.0
In [149]:
cov(X, dims=2)
Out[149]:
3×3 Matrix{Float64}:
  1.98172   -1.15707   -0.824646
 -1.15707    1.92677   -0.769697
 -0.824646  -0.769697   1.59434
In [150]:
cov(d)
Out[150]:
3×3 Matrix{Float64}:
  2.5  -1.5  -1.0
 -1.5   2.1  -0.6
 -1.0  -0.6   1.6
In [151]:
# 各次元の値をとる範囲
xs = 0:M

fig, ax = subplots()
cs = ax.imshow([pdf(d, [x1, x2, M - (x1 + x2)]) for x1 in xs, x2 in xs]', origin="lower")
fig.colorbar(cs)
set_options(ax, "x1", "x2", "";
                        grid=false)
Out[151]:
false
In [152]:
μ = 2.0
d = Poisson(μ)
X = rand(d, 100)
X'
Out[152]:
1×100 adjoint(::Vector{Int64}) with eltype Int64:
 1  1  7  2  2  1  1  2  2  3  1  1  4  …  1  3  1  7  2  3  2  3  2  1  0  1
In [153]:
max_val = maximum(X)
fig, ax = subplots()
ax.hist(X, bins=max_val+1, range=[0, max_val])
set_options(ax, "x", "freqency", "";
                        gridy=true)
Out[153]:
false
In [154]:
println("mean (exact) = $(mean(d)), var (exact) = $(var(d))")
mean (exact) = 2.0, var (exact) = 2.0
In [155]:
println("mean (approx) = $(mean(X)), var (approx) = $(var(X))")
mean (approx) = 2.05, var (approx) = 2.0277777777777777
In [156]:
# 表示範囲は0から25までとする
xs = 0:25

# 平均パラメータのリスト
μs = [0.5, 1.0, 4.0, 10.0]

fig, axes = subplots(2, 2, sharey=true, figsize=(8,6))
for (i, ax) in enumerate(axes)
    μ = μs[i]
    d = Poisson(μ)
    ax.bar(xs, pdf.(d, xs))
    set_options(ax, "x", "probability", "μ=$(μ)";
                            gridy=true)
end
tight_layout()
In [157]:
# 負の二項分布の作成
r = 10
μ = 0.3
d = NegativeBinomial(r, μ)

#値のサンプリング
X = rand(d, 100)
X'
Out[157]:
1×100 adjoint(::Vector{Int64}) with eltype Int64:
 13  33  17  34  11  29  20  36  18  …  14  20  17  24  22  24  13  24  18
In [158]:
println("mean (exact) = $(mean(d)), var (exact) = $(var(d))")
mean (exact) = 23.333333333333336, var (exact) = 77.77777777777779
In [159]:
println("maen(approx) = $(mean(X)), var (exact) = $(var(X))")
maen(approx) = 22.68, var (exact) = 67.81575757575757
In [160]:
# 表示範囲は0から60までとする
xs = 0:60

# パラメータのリスト
rs = [3, 5, 10]
μs = [0.3, 0.5, 0.7]

fig, axes = subplots(3, 3, sharey=true, figsize=(12,12))
for (i, r) in enumerate(rs)
    for (j, μ) in enumerate(μs)
        d = NegativeBinomial(r, μ)
        axes[i, j].bar(xs, pdf.(d, xs))
        
        # 平均と分散の計算、表示は小数3桁に丸める
        m = round(mean(d), digits=3)
        v = round(var(d), digits=3)
        
        axes[i, j].text(20, 0.3, "mean=$(m), var=$(v)")
        set_options(axes[i, j], "x", "probability", "r=$(r), μ=$(μ)";
                                gridy=true)
    end
end
tight_layout()
In [161]:
# パラメータ
a = 0
b = 1

# 一様分布の作成
d = Uniform(a, b)

#  サンプリング
X = rand(d, 1000)
X'
Out[161]:
1×1000 adjoint(::Vector{Float64}) with eltype Float64:
 0.0236739  0.252815  0.664046  0.892596  …  0.338695  0.00607204  0.631941
In [162]:
fig, ax = subplots()
ax.hist(X)
set_options(ax, "x", "frequency", "histogram";
                        gridy=true)
Out[162]:
false
In [163]:
# 平均値パラメータ
μ = 0.0

# 標準偏差パラメータ
σ = 1.0

# 正規分布の作成
d = Normal(μ, σ)

# サンプリング
X = rand(d, 10000)
X'
Out[163]:
1×10000 adjoint(::Vector{Float64}) with eltype Float64:
 1.42876  -1.34033  0.559439  -1.42868  …  0.162162  -0.526015  2.3709
In [164]:
fig, ax = subplots()
ax.hist(X)
set_options(ax, "x", "frequency", "histogram";
                        gridy=true)
Out[164]:
false
In [165]:
println("mean (exact) = $(mean(d)), std (exact) = $(std(d))")
println("mean (approx) = $(mean(X)), std (approx) = $(std(X))")
mean (exact) = 0.0, std (exact) = 1.0
mean (approx) = -0.009697113716219372, std (approx) = 0.9908933036261273
In [166]:
# 表示の範囲は−4から4までとする
xs = range(-4, 4, length=100)

# 平均パラメータのリスト
μs = [-1.0, 0.0, 1.0]

# 標準偏差パラメータリスト
σs = [0.3, 1.0, 1.5]

fig, axes = subplots(length(μs), length(σs), sharey = true, figsize=(8, 8))
for (i, μ) in enumerate(μs)
    for (j, σ) in enumerate(σs)
        d = Normal(μ, σ)
        axes[i, j].plot(xs, pdf.(d, xs))
        set_options(axes[i, j], "x", "density", "μ = $(μ), σ = $(σ)")
    end
end
tight_layout()
In [167]:
μ = 0.0
σ = 0.2
d = Normal(μ, σ)
pdf(d, 0.1)
Out[167]:
1.7603266338214976
In [168]:
xs = range(-1, 1, length=100)
fig, axes = subplots(2, 1, figsize=(4,8))

# 正規分布の確率密度関数をプロット
axes[1].plot(xs, pdf.(d, xs))
set_options(axes[1], "x", "density", "")

# 累積分布関数をプロット
axes[2].plot(xs, cdf.(d, xs))
set_options(axes[2], "x", "cumulation", "")

tight_layout()
In [169]:
cdf(d, 0.2) - cdf(d, 0.0)
Out[169]:
0.34134474606854304
In [170]:
X = rand(d, 10000)

# 0.0から0.2にに入ったサンプルの割合を求める
mean(0.0 .< X .< 0.2)
Out[170]:
0.3421
In [171]:
μ = 0.0
σ = 1.0
d = Normal(μ, σ)
Out[171]:
Normal{Float64}(μ=0.0, σ=1.0)
In [172]:
X_obs = [0.1, -0.1, 0.2, 0.5]
pdf.(d, X_obs)
Out[172]:
4-element Vector{Float64}:
 0.3969525474770118
 0.3969525474770118
 0.3910426939754559
 0.3520653267642995
In [173]:
prod(pdf.(d, X_obs))
Out[173]:
0.021693249867975634
In [174]:
lp = sum(logpdf.(d, X_obs))
println("logpdf = $(lp)")

# 元のpdfに戻す
println("pdf = $(exp(lp))")
logpdf = -3.830754132818691
pdf = 0.021693249867975627
In [175]:
μ = [0.0, 0.0]
Σ = [1.0 0.0
        0.0 1.0]
d = MvNormal(μ, Σ)
X = rand(d, 1000)
X
Out[175]:
2×1000 Matrix{Float64}:
 -1.57857  -0.877997  0.678171  -1.71418   …   0.293371  0.277071    2.55895
  1.38254   1.35867   0.315322   0.018942     -0.87999   0.0140066  -0.510717
In [176]:
fig, ax = subplots(figsize=(4,4))
ax.scatter(X[1,:], X[2,:], alpha=0.5)
set_options(ax, "x1", "x2", "scatter")
Out[176]:
false
In [177]:
x1s = range(-3, 3, length=100)
x2s = range(-3, 3, length=100)

fig, ax = subplots(figsize=(4,4))
cs = ax.contour(x1s, x2s, [pdf(d, [x1, x2]) for x1 in x1s, x2 in x2s]')
ax.clabel(cs, inline=true)
set_options(ax, "x1", "x2", "density")
Out[177]:
false
In [178]:
μ = [0.0, 0.0]
Σ = [1.0  0.5
        0.5 1.0]
d = MvNormal(μ, Σ)
Out[178]:
FullNormal(
dim: 2
μ: [0.0, 0.0]
Σ: [1.0 0.5; 0.5 1.0]
)
In [179]:
X = rand(d, 1000)
fig, ax = subplots(figsize=(4,4))
ax.scatter(X[1,:], X[2,:], alpha=0.5)
set_options(ax, "x1", "x2", "scatter")
Out[179]:
false
In [180]:
x1s = range(-3, 3, length=100)
x2s = range(-3, 3, length=100)

fig, ax = subplots(figsize=(4,4))
cs = ax.contour(x1s, x2s, [pdf(d, [x1, x2]) for x1 in x1s, x2 in x2s]')
ax.clabel(cs, inline=true)
set_options(ax, "x1", "x2", "density")
Out[180]:
false
In [181]:
μ = [0.0, 0.0]
Σ = [1.0  -0.5
        -0.5 1.0]
d = MvNormal(μ, Σ)
Out[181]:
FullNormal(
dim: 2
μ: [0.0, 0.0]
Σ: [1.0 -0.5; -0.5 1.0]
)
In [182]:
X = rand(d, 1000)
fig, ax = subplots(figsize=(4,4))
ax.scatter(X[1,:], X[2,:], alpha=0.5)
set_options(ax, "x1", "x2", "scatter")
Out[182]:
false
In [183]:
x1s = range(-3, 3, length=100)
x2s = range(-3, 3, length=100)

fig, ax = subplots(figsize=(4,4))
cs = ax.contour(x1s, x2s, [pdf(d, [x1, x2]) for x1 in x1s, x2 in x2s]')
ax.clabel(cs, inline=true)
set_options(ax, "x1", "x2", "density")
Out[183]:
false
In [184]:
μ = [0.0, 0.0]
Σ = [1.5  0.25
        0.25 0.5]
d = MvNormal(μ, Σ)
Out[184]:
FullNormal(
dim: 2
μ: [0.0, 0.0]
Σ: [1.5 0.25; 0.25 0.5]
)
In [185]:
function approx_integration(x_range, f)
    # 幅
    Δ = x_range[2] - x_range[1]
    
    # 近似された面積と幅を返す
    sum([f(x) * Δ for x in x_range]), Δ
end
Out[185]:
approx_integration (generic function with 1 method)
In [186]:
# 同時分布
fig, axes = subplots(2, 2, figsize=(8,8))

cs = axes[1,1].contour(x1s, x2s, [pdf(d, [x1, x2]) for x1 in x1s, x2 in x2s]')
axes[1,1].clabel(cs, inline=true)
set_options(axes[1,1], "x1", "x2", "p(x1, x2)")

# x1の周辺分布
x_range = range(-3, 3, length=100)
p1_marginal(x1) = approx_integration(x_range, x2 -> pdf(d, [x1, x2]))[1]
axes[2,1].plot(x1s, p1_marginal.(x1s))
axes[2,1].set_xlim([minimum(x1s), maximum(x1s)])
set_options(axes[2,1], "x1", "density", "p(x1)")

# x2の周辺分布
x_range = range(-3, 3, length=100)
p2_marginal(x2) = approx_integration(x_range, x1 -> pdf(d, [x1, x2]))[1]
axes[1,2].plot(p2_marginal.(x2s), x2s)
axes[1,2].set_ylim([minimum(x2s), maximum(x2s)])
set_options(axes[1,2], "density", "x2","p(x2)")

axes[2,2].axis("off")

tight_layout()
In [187]:
# 同時分布
fig, axes = subplots(2, 2, figsize=(8,8))

cs = axes[1,1].contour(x1s, x2s, [pdf(d, [x1, x2]) for x1 in x1s, x2 in x2s]')
axes[1,1].clabel(cs, inline=true)
set_options(axes[1,1], "x1", "x2", "p(x1, x2)")

# x1の条件付き分布
x2 = 1.0
p1_conditional(x1) = pdf(d, [x1, x2]) / p2_marginal(x2)
axes[1,1].plot([minimum(x1s), maximum(x1s)], [x2, x2], "r--")
axes[2,1].plot(x1s, p1_conditional.(x1s), "r", label="p(x1|x2=$(x2))")

x2 = 0.0
p1_conditional(x1) = pdf(d, [x1, x2]) / p2_marginal(x2)
axes[1,1].plot([minimum(x1s), maximum(x1s)], [x2, x2], "g--")
axes[2,1].plot(x1s, p1_conditional.(x1s), "g", label="p(x1|x2=$(x2))")

axes[2,1].set_xlim([minimum(x1s), maximum(x1s)],)
set_options(axes[2,1], "x1", "density", "p(x1|x2)")

# x2の条件付き分布
x1 = 2.0
p2_conditional(x2) = pdf(d, [x1, x2]) / p1_marginal(x1)
axes[1,1].plot([x1,x1], [minimum(x2s), maximum(x2s)], "r--")
axes[1,2].plot(p2_conditional.(x2s), x2s, "r", label="p(x2|x1=$(x1))")

x1 = -2.0
p2_conditional(x2) = pdf(d, [x1, x2]) / p1_marginal(x1)
axes[1,1].plot([x1,x1], [minimum(x2s), maximum(x2s)], "g--")
axes[1,2].plot(p2_conditional.(x2s), x2s, "g", label="p(x2|x1=$(x1))")

axes[1,2].set_ylim([minimum(x2s), maximum(x2s)],)
set_options(axes[1,2], "x2", "density", "p(x2|x1)")

axes[2,2].axis("off")

tight_layout()
In [188]:
# ガンマ分布を作成
α = 1.5
θ = 2.5
d = Gamma(α, θ)

# サンプルの生成
X = rand(d, 100)
X'
Out[188]:
1×100 adjoint(::Vector{Float64}) with eltype Float64:
 0.813145  2.97026  7.63536  0.301381  …  4.02372  7.50194  2.93235  4.89421
In [189]:
fig, ax = subplots()
ax.hist(X)
set_options(ax, "x", "frequency", "histogram")
Out[189]:
false
In [190]:
println("mean (exact) : $(mean(d)), $(α*θ)")
println("var (exact) : $(var(d)), $(α*θ^2)")
mean (exact) : 3.75, 3.75
var (exact) : 9.375, 9.375
In [191]:
xs = range(0, 10, length=100)
αs = [0.5, 1.0, 2.0]
θs = [0.5, 1.0, 1.5]
fig, axes = subplots(length(αs), length(θs), sharey=true, figsize=(8,8))
for (i, α) in enumerate(αs)
    for (j, θ) in enumerate(θs)
        d = Gamma(α, θ)
        μ, σ = mean(d), std(d)
        axes[i,j].plot(xs, pdf.(d, xs))
        axes[i,j].set_ylim([0, 1.5])
        set_options(axes[i,j], "x", "density", "α=$(α), θ=$(θ), \n" * 
                                                "(μ=$(round(μ, digits=3)), σ=$(round(σ, digits=3)))")
    end
end
tight_layout()
In [192]:
# ベータ分布の作成
α = 0.5
β = 0.5
d = Beta(α, β)

# サンプルの生成
X = rand(d, 100)
X'
Out[192]:
1×100 adjoint(::Vector{Float64}) with eltype Float64:
 0.610743  0.745976  0.998187  0.525525  …  0.0209475  0.315101  0.739733
In [193]:
fig, ax = subplots()
ax.hist(X)
set_options(ax, "x", "frequency", "histogram")
Out[193]:
false
In [194]:
xs = range(0, 1, length=100)

# パラメータのリスト
αs = [0.1, 1.0, 2.0]
βs = [0.1, 1.0, 2.0]

fig, axes = subplots(length(αs), length(βs), sharey=true, figsize=(8,8))
for (i, α) in enumerate(αs)
    for (j, β) in enumerate(βs)
        d = Beta(α, β)
        μ, σ = mean(d), std(d)
        axes[i,j].plot(xs, pdf.(d, xs))
        set_options(axes[i,j], "x", "density", "α=$(α), β=$(β), \n" *
                                                    "(μ=$(round(μ, digits=3)), σ=$(round(σ, digits=3)))")
    end
end
tight_layout()
In [195]:
# ディリクレ分布の作成
α = [0.75, 0.75, 0.75]
d = Dirichlet(α)

# サンプルの生成
X = rand(d, 1000)
X
Out[195]:
3×1000 Matrix{Float64}:
 0.280472  0.422914  0.108769  0.423537  …  0.176149  0.785769   0.0221444
 0.145596  0.457917  0.435515  0.160961     0.304625  0.0926957  0.977021
 0.573932  0.119168  0.455715  0.415503     0.519226  0.121535   0.000834433
In [196]:
fig, axes = subplots(1, 2, figsize=(8,4))

# 散布図による可視化
axes[1].scatter(X[1,:], X[2,:], alpha=0.25)
set_options(axes[1], "x1", "x2", "scatter")

# 確率密度関数
x1s = range(0, 1, length=100)
x2s = range(0, 1, length=100)
cs = axes[2].contour(x1s, x2s, [x1 + x2 > 1 ? 0.0 : pdf(d, [x1, x2, 1 - (x1 + x2)]) for x1 in x1s, x2 in x2s]')
axes[2].clabel(cs, inline=true)
set_options(axes[2], "x1", "x2", "density")
tight_layout()
In [197]:
mean(d)
Out[197]:
3-element Vector{Float64}:
 0.3333333333333333
 0.3333333333333333
 0.3333333333333333
In [198]:
cov(d)
Out[198]:
3×3 Matrix{Float64}:
  0.0683761  -0.034188   -0.034188
 -0.034188    0.0683761  -0.034188
 -0.034188   -0.034188    0.0683761
In [199]:
# パラメータのリスト
αs = [[0.1, 0.1, 0.1],
            [0.5, 0.5, 0.5],
            [1.0, 1.0, 1.0],
            [2.0, 2.0, 2.0],
            [5.0, 5.0, 5.0],
            [0.1, 0.1, 0.5],
            [0.1, 0.5, 0.1],
            [0.1, 0.5, 5.0],
            [1.0, 2.0, 5.0]]

xs = range(0.1, 0.99, length=100)
ys = range(0.1, 0.99, length=100)
fig, axes = subplots(3, 3, figsize=(9,9))
for (i, α) in enumerate(αs)
    print(α)
    d = Dirichlet(α)
    X = rand(d, 1000)
    axes[i].scatter(X[1,:], X[2,:], alpha=0.25)
    axes[i].set_xlim([0,1])
    axes[i].set_ylim([0,1])
    set_options(axes[i], "x1", "x2", "α=$(α)")
end
tight_layout()
[0.1, 0.1, 0.1][0.5, 0.5, 0.5][1.0, 1.0, 1.0][2.0, 2.0, 2.0][5.0, 5.0, 5.0][0.1, 0.1, 0.5]
[0.1, 0.5, 0.1][0.1, 0.5, 5.0][1.0, 2.0, 5.0]
In [200]:
μ = 0.0
σ = 1.0
d = Normal(μ, σ)
X = rand(d, 100)
Y = exp.(X)
Y'
Out[200]:
1×100 adjoint(::Vector{Float64}) with eltype Float64:
 0.856311  0.832838  3.73996  4.16099  …  2.29158  0.986543  1.66065  4.80166
In [201]:
fig, axes = subplots(1, 2, figsize=(12,4))
axes[1].hist(X)
set_options(axes[1], "x", "freqency", "Normal distribution";
                        gridy=true)
axes[2].hist(Y)
set_options(axes[2], "y", "frequency", "Log-Normal distribution";
                        gridy=true)
tight_layout()
In [202]:
d = LogNormal(μ, σ)
X = rand(d, 100)
X'
Out[202]:
1×100 adjoint(::Vector{Float64}) with eltype Float64:
 0.558211  2.9403  0.564984  0.803061  …  1.96972  1.23392  0.1873  2.36384
In [203]:
fig, ax = subplots()
ax.hist(X)
set_options(ax, "x", "freqency", "histogram";
                        gridy=true)
Out[203]:
false
In [204]:
xs = range(0, 5, length=100)
μs = [-1, 0.0, 1.0]
σs = [0.2, 1.0, 1.5]
fig, axes = subplots(length(μs), length(σs), sharey=true, figsize=(12,12))
for (i, μ) in enumerate(μs)
    for (j, σ) in enumerate(σs)
        d = LogNormal(μ, σ)
        axes[i,j].plot(xs, pdf.(d, xs))
        m = round(mean(d), digits=2)
        v = round(var(d), digits=2)
        axes[i,j].text(2, 5, "mean = $(m), var = $(v)")
        set_options(axes[i,j], "x", "density", "μ = $(μ), σ = $(σ)")
    end
end
tight_layout()
In [205]:
# パラメータを生成するための分布
μ1 = [-1.0, 1.0]
Σ1 = [0.2  0.0
            0.0 0.2]
μ2 = [1.0, -1.0]
Σ2 = [0.4 0.0
            0.0 0.4]
p = 0.3

# サンプリングする回数
N = 100

# 生成データを保存する配列
X = Array{Float64}(undef, 2, N)

# 選択された分布を示す配列(潜在変数)
S = Array{Bool}(undef, N)

for i in 1:N
    # 潜在変数のサンプル
    S[i] = rand(Bernoulli(p))
    
    # 潜在変数の値に応じて多変量正規分布のパラメータを切り替える
    (μ, Σ) = S[i] == 1 ? (μ1, Σ1) : (μ2, Σ2)
    
    # データをサンプリングする
    X[:, i] = rand(MvNormal(μ, Σ))
end

# 生成されたデータを散布図として可視化
fig, ax = subplots(figsize=(4,4))
ax.scatter(X[1,:], X[2,:], alpha=0.5)
set_options(ax, "x1", "x2", "scatter")
Out[205]:
false
In [206]:
# 混合多変量正規分布の確率密度関数
pdfgmm(x) = p*pdf(MvNormal(μ1, Σ1), x) + (1-p)*pdf(MvNormal(μ2, Σ2), x)

xs1 = range(-3, 3, length=100)
xs2 = range(-3, 3, length=100)

fig, ax = subplots(figsize=(4,4))
cs = ax.contour(xs1, xs2, [pdfgmm([x1, x2]) for x1 in xs1, x2 in xs2]')
ax.clabel(cs, inline=true)
set_options(ax, "x1", "x2", "density")
Out[206]:
false
In [207]:
# パラメータを生成するための分布
μ = [0.0, 0.0]
Σ = [0.1 0.0
        0.0 0.1]

# 出力値に付加するノイズ
σ = 1.0

# 入力値はあらかじめ与えておく
X = [-10, -5, 0, 5, 10]

# サンプリングする回数
num_samples = 3

# パラメータのサンプル
W = rand(MvNormal(μ, Σ), num_sample)

# 出力値を保存するためのリスト
Ys = []

fig, axes = subplots(1, 2, figsize=(8,3))
xs = range(-12, 12, length=100)
for n in 1:num_samples
    w1, w2 = W[:,n]
    
    # パラメータをプロット
    axes[1].scatter(w1, w2)
    
    # 生成された関数をプロット
    linear_function(x) = w1*x + w2
    axes[2].plot(xs, linear_function.(xs))
    
    # 関数からの出力値も生成する
    Y = rand.(Normal.(linear_function.(X), σ))
    push!(Ys, Y)
end
set_options(axes[1], "w1", "w2", "parameters")
set_options(axes[2], "x", "y", "functions")
tight_layout()
UndefVarError: num_sample not defined

Stacktrace:
 [1] top-level scope
   @ In[207]:16
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [208]:
fig, axes = subplots(1, num_sample, sharey=true, figsize=(12,3))
for n in 1:num_sample
    w1, w2 = W[:,n]
    Y = Ys[n]
    linear_function(x) =w1*x + w2
    axes[n].plot(xs, linear_function.(xs))
    axes[n].scatter(X, Y)
    set_options(axes[n], "x", "y", "w1 = $(round(w1, digits=3)), w2 = $(round(w2, digits=3))")
end
tight_layout()
UndefVarError: num_sample not defined

Stacktrace:
 [1] top-level scope
   @ In[208]:1
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [209]:
sig(x) = 1 / (1 + exp(-x))
Out[209]:
sig (generic function with 1 method)
In [210]:
# パラメータを生成するための分布
μ = [0.0, 0.0]
Σ = [0.01 0.0
        0.0 0.01]

# 入力値はあらかじめ与えておく
X = [-10, -5, 0, 5, 10]

# サンプリングする回数
num_samples = 3

# パラメータのサンプル
W = rand(MvNormal(μ, Σ), num_sample)

# 出力値を保存するためのリスト
Ys = []

fig, axes = subplots(1, 2, figsize=(8,3))
xs = range(-12, 12, length=100)
for n in 1:num_samples
    w1, w2 = W[:,n]
    
    # パラメータをプロット
    axes[1].scatter(w1, w2)
    
    # 生成された関数をプロット
    nonlinear_function(x) = sig(w1*x + w2)
    axes[2].plot(xs, nonlinear_function.(xs))
    
    # 関数からの出力値も生成する
    Y = rand.(Bernoulli.(nonlinear_function.(X)))
    push!(Ys, Y)
end
set_options(axes[1], "w1", "w2", "parameters")
set_options(axes[2], "x", "y", "functions")
tight_layout()
UndefVarError: num_sample not defined

Stacktrace:
 [1] top-level scope
   @ In[210]:13
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [211]:
fig, axes = subplots(1, num_sample, sharey=true, figsize=(12,3))
xs = range(-12, 12, length=100)
for n in 1:num_samples
    w1, w2 = W[:,n]
    Y = Ys[n]
    nonlinear_function(x) = sig(w1*x + w2)
    axes[n].plot(xs, nonlinear_function.(xs))
    axes[n].scatter(X, Y)
    set_options(axes[n], "x", "y", "w1=$(round(w1, digits=3)), w2=$(round(w2, digits=3))")
end
tight_layout()
UndefVarError: num_sample not defined

Stacktrace:
 [1] top-level scope
   @ In[211]:1
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [212]:
# パラメータを生成するための分布
μ = [0.0, 0.0]
Σ = [0.01 0.0
        0.0 0.01]

# 入力値はあらかじめ与えておく
X = [-10, -5, 0, 5, 10]

# サンプリングする回数
num_samples = 3

# パラメータのサンプル
W = rand(MvNormal(μ, Σ), num_sample)

# 出力値を保存するためのリスト
Ys = []

fig, axes = subplots(1, 2, figsize=(8,3))
xs = range(-12, 12, length=100)
for n in 1:num_samples
    w1, w2 = W[:,n]
    
    # パラメータをプロット
    axes[1].scatter(w1, w2)
    
    # 生成された関数をプロット
    nonlinear_function(x) = exp(w1*x + w2)
    axes[2].plot(xs, nonlinear_function.(xs))
    
    # 関数からの出力値も生成する
    Y = rand.(Poisson.(nonlinear_function.(X)))
    push!(Ys, Y)
end
set_options(axes[1], "w1", "w2", "parameters")
set_options(axes[2], "x", "y", "functions")
tight_layout()
UndefVarError: num_sample not defined

Stacktrace:
 [1] top-level scope
   @ In[212]:13
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [213]:
fig, axes = subplots(1, num_sample, sharey=true, figsize=(12,3))
xs = range(-12, 12, length=100)
for n in 1:num_samples
    w1, w2 = W[:,n]
    Y = Ys[n]
    nonlinear_function(x) = exp(w1*x + w2)
    axes[n].plot(xs, nonlinear_function.(xs))
    axes[n].scatter(X, Y)
    set_options(axes[n], "x", "y", "w1=$(round(w1, digits=3)), w2=$(round(w2, digits=3))")
end
tight_layout()
UndefVarError: num_sample not defined

Stacktrace:
 [1] top-level scope
   @ In[213]:1
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [214]:
# パッケージの読み込み
using Distributions, PyPlot, LinearAlgebra
WARNING: using LinearAlgebra.I in module Main conflicts with an existing identifier.
In [215]:
function generate(N)
    μ = rand(Uniform(0, 1))
    X = rand(Bernoulli(μ), N)
    μ, X
end

# 5回コイン投げを行う
generate(5)
Out[215]:
(0.14466576460977287, Bool[0, 0, 0, 1, 0])
In [216]:
# 1を表、0を裏とする
side(x) = x == 1 ? "表" : "裏"

for i in 1:10
    μ, X = generate(5)
    println("コイン$(i)、表が出る確率 μ = $(μ)、出目 X = $(side.(X))")
end
コイン1、表が出る確率 μ = 0.2115239161862532、出目 X = ["表", "表", "表", "裏", "裏"]
コイン2、表が出る確率 μ = 0.19157334256084235、出目 X = ["表", "表", "裏", "表", "表"]
コイン3、表が出る確率 μ = 0.6340337570025956、出目 X = ["表", "表", "表", "裏", "表"]
コイン4、表が出る確率 μ = 0.4673419080511929、出目 X = ["裏", "表", "表", "表", "裏"]
コイン5、表が出る確率 μ = 0.7624851169457685、出目 X = ["表", "裏", "表", "裏", "裏"]
コイン6、表が出る確率 μ = 0.7451980753784295、出目 X = ["表", "裏", "表", "裏", "表"]
コイン7、表が出る確率 μ = 0.26375329479765175、出目 X = ["表", "裏", "表", "裏", "表"]
コイン8、表が出る確率 μ = 0.8119528312066475、出目 X = ["裏", "表", "表", "表", "表"]
コイン9、表が出る確率 μ = 0.6140591793610704、出目 X = ["表", "裏", "裏", "表", "表"]
コイン10、表が出る確率 μ = 0.33996595666852136、出目 X = ["裏", "表", "裏", "裏", "裏"]
In [217]:
μs = range(0, 1, length=100)
fig, ax = subplots()
ax.plot(μs, pdf.(Uniform(0, 1), μs))
set_options(ax, "μ", "density", "Uniform distribution")
ax.set_xlim([0, 1])
ax.set_ylim([0, 1.1])
Out[217]:
(0.0, 1.1)
In [218]:
# 「裏、裏、裏、表、表」をデータとして取得
X_obs1 = [0,0,0,1,1]
Out[218]:
5-element Vector{Int64}:
 0
 0
 0
 1
 1
In [219]:
maxiter = 1_000_000
μ_posteriori1 = []
for i in 1:maxiter
    # パラメータおよびデータの生成
    μ, X = generate(length(X_obs1))
    
    #  X内の1の合計が観測と一致していれば、この時のパラメータを受容
    sum(X) == sum(X_obs1) && push!(μ_posteriori1, μ)
end

# 受容率の計算
acceptance_rate = length(μ_posteriori1) / maxiter
println("acceptance rate = $(acceptance_rate)")

μ_posteriori'
acceptance rate = 0.166686
UndefVarError: μ_posteriori not defined

Stacktrace:
 [1] top-level scope
   @ In[219]:15
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [220]:
fig, ax = subplots()
ax.hist(μ_posteriori1)
ax.set_xlim([0,1])
set_options(ax, "μ", "frequency", "histogram";
                        gridy=true)
Out[220]:
false
In [221]:
X_obs2 = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1]
Out[221]:
20-element Vector{Int64}:
 0
 0
 0
 0
 0
 0
 0
 0
 0
 0
 0
 0
 1
 1
 1
 1
 1
 1
 1
 1
In [222]:
maxiter = 1_000_000
μ_posteriori2 = []
for i in 1:maxiter
    # パラメータおよびデータの生成
    μ, X = generate(length(X_obs2))
    
    #  X内の1の合計が観測と一致していれば、この時のパラメータを受容
    sum(X) == sum(X_obs2) && push!(μ_posteriori2, μ)
end

# 受容率の計算
acceptance_rate = length(μ_posteriori2) / maxiter
println("acceptance rate = $(acceptance_rate)")

μ_posteriori'
acceptance rate = 0.047841
UndefVarError: μ_posteriori not defined

Stacktrace:
 [1] top-level scope
   @ In[222]:15
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [223]:
fig, ax = subplots()
ax.hist(μ_posteriori2)
ax.set_xlim([0,1])
set_options(ax, "μ", "frequency", "histogram";
                        gridy=true)
Out[223]:
false
In [224]:
pred1 = mean(rand.(Bernoulli.(μ_posteriori1)))
pred2 = mean(rand.(Bernoulli.(μ_posteriori2)))

println("$(pred1), $(pred2)")
0.4280563454639262, 0.40691039066909135
In [225]:
function generate2(N)
    # μの事前分布を修正
    μ = rand(Uniform(0, 0.5))
    X = rand(Bernoulli(μ), N)
    μ, X
end
generate2(5)
Out[225]:
(0.33987893554338955, Bool[0, 0, 0, 0, 1])
In [226]:
for i in 1:10
    μ, X = generate2(5)
    println("コイン $(i)、表が出る確率 μ = $(μ)、出目 X = $(side.(X))")
end
コイン 1、表が出る確率 μ = 0.371291143488243、出目 X = ["表", "裏", "裏", "表", "裏"]
コイン 2、表が出る確率 μ = 0.11257772410611666、出目 X = ["裏", "裏", "裏", "表", "表"]
コイン 3、表が出る確率 μ = 0.4383712405951564、出目 X = ["表", "裏", "表", "表", "裏"]
コイン 4、表が出る確率 μ = 0.23543633208856873、出目 X = ["裏", "裏", "裏", "裏", "裏"]
コイン 5、表が出る確率 μ = 0.0069485910271677165、出目 X = ["裏", "裏", "裏", "裏", "裏"]
コイン 6、表が出る確率 μ = 0.05611836072281462、出目 X = ["裏", "裏", "裏", "裏", "裏"]
コイン 7、表が出る確率 μ = 0.4629175852579448、出目 X = ["表", "表", "裏", "表", "裏"]
コイン 8、表が出る確率 μ = 0.35259901883072686、出目 X = ["裏", "表", "表", "表", "裏"]
コイン 9、表が出る確率 μ = 0.03320618117563745、出目 X = ["裏", "裏", "裏", "裏", "裏"]
コイン 10、表が出る確率 μ = 0.15212636716313133、出目 X = ["裏", "裏", "裏", "裏", "裏"]
In [227]:
maxiter = 1_000_000
μ_posteriori3 = []
for i in 1:maxiter
    # パラメータおよびデータの生成
    μ, X = generate2(length(X_obs1))
    
    #  X内の1の合計が観測と一致していれば、この時のパラメータを受容
    sum(X) == sum(X_obs1) && push!(μ_posteriori3, μ)
end

# 受容率の計算
acceptance_rate = length(μ_posteriori3) / maxiter
println("acceptance rate = $(acceptance_rate)")

μ_posteriori'
acceptance rate = 0.219677
UndefVarError: μ_posteriori not defined

Stacktrace:
 [1] top-level scope
   @ In[227]:15
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116
In [228]:
fig, ax = subplots()
ax.hist(μ_posteriori3)
ax.set_xlim([0,1])
set_options(ax, "μ", "frequency", "histogram";
                        gridy=true)
Out[228]:
false
In [229]:
# 同時分布p(X, μ)の確率密度関数
p_joint(X, μ) = prod(pdf.(Bernoulli(μ), X)) * pdf(Uniform(0,1), μ)

#  数値積分
function approx_integration(μ_range, p)
    Δ = μ_range[2] - μ_range[1]
    X -> sum([p(X, μ) * Δ for μ in μ_range]), Δ
end

# μの積分範囲
μ_range = range(0, 1, length=100)

# 数値積分の実行
p_marginal, Δ = approx_integration(μ_range, p_joint)

#  データ(2種類)
X_obs1 = [0,0,0,1,1]
X_obs2 = [0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1]

# それぞれの周辺尤度の近似計算
println("$(p_marginal(X_obs1)), $(p_marginal(X_obs2))")
0.016666666493163274, 3.7801895387034807e-7
In [230]:
# パラメータの可視化範囲
μs = range(0, 1, length=100)

fig, axes = subplots(1, 2, sharey=true, figsize=(10,4))
for (i, X_obs) in enumerate([X_obs1, X_obs2])
    posterior(μ) = p_joint(X_obs, μ) / p_marginal(X_obs)
    axes[i].plot(μs, posterior.(μs))
    set_options(axes[i], "μ", "density", "approximate posterior (X_obs$(i))")
end
In [231]:
# 積分の中身の式
posterior1(μ) = p_joint(X_obs1, μ) / p_marginal(X_obs1)
posterior2(μ) = p_joint(X_obs2, μ) / p_marginal(X_obs2)
p_inner1(x, μ) = pdf.(Bernoulli(μ), x) * posterior1(μ)
p_inner2(x, μ) = pdf.(Bernoulli(μ), x) * posterior2(μ)

# パラメータμに関する積分
μ_range = range(0, 1, length=100)
pred1, Δ1 = approx_integration(μ_range, p_inner1)
pred2, Δ2 = approx_integration(μ_range, p_inner2)

# 1が出る予測確率
println("$(pred1(1)), $(pred2(1))")
0.4285714434416308, 0.40909090909090784
In [232]:
fig, axes = subplots(1, 2, sharey=true, figsize=(10,4))
μs = range(0, 1, length=100)
for (i, X_obs) in enumerate([X_obs1, X_obs2])
    # 厳密な事後分布はベータ分布
    α = 1.0 + sum(X_obs)
    β = 1.0 + length(X_obs) - sum(X_obs)
    d = Beta(α, β)
    
    # 事後分布を可視化
    axes[i].plot(μs, pdf.(d, μs))
    set_options(axes[i], "μ", "density", "exact posterior (X_obs$(i))")
end
In [233]:
function prediction(X_obs)
    α = 1.0 + sum(X_obs)
    β = 1.0 + length(X_obs) - sum(X_obs)
    α / (α + β)
end

println("$(prediction(X_obs1)), $(prediction(X_obs2))")
0.42857142857142855, 0.4090909090909091
In [234]:
function generate_linear(X, σ, μ1, μ2, σ1, σ2)
    w1 = rand(Normal(μ1, σ1))
    w2 = rand(Normal(μ2, σ2))
    linear_function(x) = w1*x + w2
    Y = rand.(Normal.(linear_function.(X), σ))
    Y, linear_function, w1, w2
end
Out[234]:
generate_linear (generic function with 1 method)
In [235]:
# あらかじめ与えられるパラメータと入力集合X
σ = 1.0
μ1 = 0.0
μ2 = 0.0
σ1 = 10.0
σ2 = 10.0
X = [-1.0, -0.5, 0, 0.5, 1.0]

# 可視化する範囲
xs = range(-2, 2, length=100)

fig, axes = subplots(2, 3, sharey=true, figsize=(12,6))
for ax in axes
    # 関数f、出力Yの生成
    Y, linear_function, w1, w2 = generate_linear(X, σ, μ1, μ2, σ1, σ2)
    
    # 生成された直線とYのプロット
    ax.plot(xs, linear_function.(xs), label="simulated function")
    ax.scatter(X, Y, label="simulated data")
    
    set_options(ax, "x", "y", "data (N = $(length(X)))", legend=true)
end
tight_layout()
In [236]:
# 平均パラメータのリスト
μ1_list = [-20.0, 0.0, 20.0]
μ2_list = [-20.0, 0.0, 20.0]

# 標準偏差パラメータは固定
σ1 = 10.0
σ2 = 10.0

fig, axes = subplots(length(μ1_list), length(μ2_list), sharey=true, figsize=(12, 12))
for (i, μ1) in enumerate(μ1_list)
    for(j, μ2) in enumerate(μ2_list)
        # 関数を複数サンプル
        fs = [generate_linear(X, σ, μ1, μ2, σ1, σ2)[2] for _ in 1:100]
        
        # 生成された直線群をプロット
        for linear_function in fs
            axes[i,j].plot(xs, linear_function.(xs), "g", alpha=0.1)
        end
        
        axes[i,j].set_xlim(extrema(xs))
        set_options(axes[i,j], "x", "y", "μ1=$(μ1), μ2=$(μ2)")
    end
end
tight_layout()
In [237]:
# 標準偏差パラメータのリスト
σ1_list = [1.0, 10.0, 20.0]
σ2_list = [1.0, 10.0, 20.0]

# 平均パラメータは固定
μ1 = 0.0
μ2 = 0.0

fig, axes = subplots(length(σ1_list), length(σ2_list), sharey=true, figsize=(12, 12))
for (i, σ1) in enumerate(σ1_list)
    for(j, σ2) in enumerate(σ2_list)
        # 関数を複数サンプル
        fs = [generate_linear(X, σ, μ1, μ2, σ1, σ2)[2] for _ in 1:100]
        
        # 生成された直線群をプロット
        for linear_function in fs
            axes[i,j].plot(xs, linear_function.(xs), "g", alpha=0.1)
        end
        
        axes[i,j].set_xlim(extrema(xs))
        set_options(axes[i,j], "x", "y", "σ1=$(σ1), σ2=$(σ2)")
    end
end
tight_layout()
In [238]:
# 入力データ
X_obs = [-2, 1, 5]

# 出力データ
Y_obs = [-2.2, -1.0, 1.5]

# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
set_options(ax, "x", "y", "data = (N = $(length(X_obs)))")
Out[238]:
false
In [239]:
# 同時分布
p_joint(X, Y, w1, w2) = prod(pdf.(Normal.(w1.*X .+ w2, σ), Y)) * pdf(Normal(μ1, σ1), w1) * pdf(Normal(μ2, σ2), w2)

# 数値積分
function approx_integration_2D(w_range, p)
    Δ = w_range[2] - w_range[1]
    (X, Y) -> sum([p(X, Y, w1, w2) * Δ^2 for w1 in w_range, w2 in w_range])
end

# wの積分範囲
w_range = range(-3, 3, length=100)

# 数値積分の実行
p_marginal = approx_integration_2D(w_range, p_joint)
p_marginal(X_obs, Y_obs)
Out[239]:
6.924264340150274e-5
In [240]:
# 事後分布の計算
w_posterior = [p_joint(X_obs, Y_obs, w1, w2) for w1 in w_range, w2 in w_range] ./ p_marginal(X_obs, Y_obs)

fig, axes = subplots(1, 2, figsize=(8,4))

# 等高線図
cs = axes[1].contour(w_range, w_range, w_posterior', cmap="jet")
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "posterior density (contour)")

# カラーメッシュ
xgrid = repeat(w_range', length(w_range), 1)
ygrid = repeat(w_range, 1, length(w_range))
axes[2].pcolormesh(xgrid, ygrid, w_posterior', cmap="jet", shading="auto")
set_options(axes[2], "w1", "w2", "posterior density (colored)")

tight_layout()
In [241]:
function approx_predictive(w_posterior, w_range, p)
    Δ = w_range[2] - w_range[1]
    (x, y) -> sum([p(x, y, w1, w2) * w_posterior[i,j] * Δ^2 for (i, w1) in enumerate(w_range), (j, w2) in enumerate(w_range)])
end
p_likelihood(xp, yp, w1, w2) = pdf(Normal(w1*xp + w2, σ), yp)
p_predictive = approx_predictive(w_posterior, w_range, p_likelihood)
Out[241]:
#80 (generic function with 1 method)
In [242]:
xp = 4.0
fig, ax = subplots()
ys = range(-5, 5, length=100)
ax.plot(ys, p_predictive.(xp, ys))
set_options(ax, "y_p", "density", "predictive distiribution p(y_p|x_p=$(xp), X=X_obs, Y=Y_obs)")
Out[242]:
false
In [243]:
# 描画範囲
xs = range(-10, 10, length=100)
ys = range(-5, 5, length=100)

# 密度の計算
density_y = p_predictive.(xs, ys')

fig, axes = subplots(1, 2, sharey=true, figsize=(8,4))

# 等高線図
cs = axes[1].contour(xs, ys, density_y', cmap="jet")
axes[1].clabel(cs, inline=true)
axes[1].scatter(X_obs, Y_obs)
set_options(axes[1], "x", "y", "predictive distribution (contour)")

# カラーメッシュ
xgrid = repeat(xs', length(ys), 1)
ygrid = repeat(ys, 1, length(xs))
axes[2].pcolormesh(xgrid, ygrid, density_y', cmap="jet", shading="auto")
axes[2].plot(X_obs, Y_obs, "ko", label="data")
set_options(axes[2], "x", "y", "predictive distribution (colored)")

tight_layout()
In [244]:
# 切片用の擬似データ
X_extended = hcat(X_obs, ones(size(X_obs)))'

# 事前分布のパラメータ
Σ = [σ1^2 0
        0 σ2^2]
μ = [μ1, μ2]

# 事後分布のパラメータ
Σ_hat = inv(σ^(-2) * X_extended*X_extended' + inv(Σ))
μ_hat = Σ_hat*(σ^(-2) * X_extended*Y_obs + inv(Σ)*μ)

# 予測分布のパラメータ
μp(xp) = (μ_hat' * xp)[1]
σp(xp) = sqrt(σ^2 + (xp' * Σ_hat * xp)[1])

# 予測分布の確率密度関数
p_predictive_exact(xp, yp) = pdf(Normal(μp(xp), σp(xp)), yp)

# 描画範囲
xs = range(-10, 10, length=100)
ys = range(-5, 5, length=100)
density_y = [p_predictive_exact([x, 1.0], y) for x in xs, y in ys]

fig, axes = subplots(1, 2, sharey=true, figsize=(8,4))

# 等高線図
cs = axes[1].contour(xs, ys, density_y', cmap="jet")
axes[1].clabel(cs, inline=true)
axes[1].scatter(X_obs, Y_obs)
set_options(axes[1], "x", "y", "predictive distiribution (contour)")

# カラーメッシュ
xgrid = repeat(xs', length(ys), 1)
ygrid = repeat(ys, 1, length(xs))
axes[2].pcolormesh(xgrid, ygrid, density_y', cmap="jet", shading="auto")
axes[2].plot(X_obs, Y_obs, "ko", label="data")
set_options(axes[2], "x", "y", "predictive distribution (colored)")

tight_layout()
In [245]:
# シグモイド関数を定義
sig(x) = 1/(1+exp(-x))

# パラメータ、関数、出力集合Yを生成
function generate_logistic(X, μ1, μ2, σ1, σ2)
    w1 = rand(Normal(μ1, σ1))
    w2 = rand(Normal(μ2, σ2))
    nonlinear_function(x) = sig(w1*x + w2)
    Y = rand.(Bernoulli.(nonlinear_function.(X)))
    Y, nonlinear_function, w1, w2
end
Out[245]:
generate_logistic (generic function with 1 method)
In [246]:
# あらかじめ与えられたパラメータと入力X
μ1 = 0
μ2 = 0
σ1 = 10.0
σ2 = 10.0
X = [-1.0, -0.5, 0, 0.5, 1.0]

# 可視化する範囲
xs = range(-2, 2, length=100)

fig, axes = subplots(2, 3, sharey=true, figsize=(12,6))
for ax in axes
    # 関数f、出力Yの生成
    Y, nonlinear_function, w1, w2 = generate_logistic(X, μ1, μ2, σ1, σ2)
    
    # 生成された直線とYのプロット
    ax.plot(xs, nonlinear_function.(xs), label="simulated function")
    ax.scatter(X, Y, label="simulated data")
    
    set_options(ax, "x", "y", "data (N = $(length(X)))", legend=true)
end
tight_layout()
In [247]:
# 入力データセット
X_obs = [-2, 1, 2]

# 出力データセット
Y_obs = Bool.([0, 1, 1])

# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
ax.set_xlim([-3, 3])
set_options(ax, "x", "y", "data (N = $(length(X_obs)))")
Out[247]:
false
In [248]:
# 最大サンプリング数
maxiter = 10_000

# パラメータ保存用
param_posterior = Vector{Tuple{Float64, Float64}}()

for i in 1:maxiter
    # 関数f、出力Yの生成
    Y, nonlinear_function, w1, w2 = generate_logistic(X_obs, μ1, μ2, σ1, σ2)
    
    # 観測データと一致していれば、そのときのパラメータwを保存
    Y == Y_obs && push!(param_posterior, (w1, w2))
end

# サンプル受容率
acceptance_rate = length(param_posterior) / maxiter
println("acceptance_rate = $(acceptance_rate)")
acceptance_rate = 0.2979
In [249]:
# パラメータ抽出用の関数
unzip(a) = map(x -> getfield.(a, x), fieldnames(eltype(a)))

# 事前分布からのサンプル(10,000個)
param_prior = [generate_logistic(X, μ1, μ2, σ1, σ2)[3:4] for i in 1:10_000]
w1_prior, w2_prior = unzip(param_prior)

# 事後分布からのサンプル
w1_posterior, w2_posterior = unzip(param_posterior)

fig, axes = subplots(1, 2, sharex=true, sharey=true, figsize=(8,4))

# 事前分布
axes[1].scatter(w1_prior, w2_prior, alpha=0.01)
set_options(axes[1], "w1", "w2", "sample from prior")

# 事後分布
axes[2].scatter(w1_posterior, w2_posterior, alpha=0.01)
set_options(axes[2], "w1", "w2", "sample from posterior")

tight_layout()
In [250]:
# 関数を可視化する範囲
xs = range(-3, 3, length=100)

# サンプリングで得られたパラメータ全体のプロット
fig, ax = subplots()
ax.scatter(w1_posterior, w2_posterior, alpha=0.1)
set_options(ax, "w1", "w2", "sample from posterior")

fig, axes = subplots(2, 3, figsize=(12,8))
for i in eachindex(axes)
    # 関数を可視化するためのwを1つ適当に選択
    j = round(Int, length(param_posterior)*rand()) + 1
    w1, w2 = param_posterior[j]
    
    # 選択されたw
    ax.scatter(w1, w2, color="r")
    ax.text(w1, w2, i)
    
    # 対応する関数のプロット
    nonlinear_function(x) = sig(w1*x + w2)
    axes[i].plot(xs, nonlinear_function.(xs), "r")
    
    # 観測データのプロット
    axes[i].scatter(X_obs, Y_obs)
    
    axes[i].set_ylim([-0.1, 1.1])
    set_options(axes[i], "x", "y (prob.)", "($(i)) w1=$(round(w1, digits=3)), w2=$(round(w2, digits=3))")
end
tight_layout()
In [251]:
fig, axes = subplots(1, 2, sharey=true, figsize=(12,4))

fs = []
for (i, param) in enumerate(param_posterior)
    # 1サンプル分のパラメータ
    w1, w2 = param
    
    # 1サンプル分の予測関数
    nonlinear_function(x) = sig(w1*x + w2) 
    axes[1].plot(xs, nonlinear_function.(xs), "g", alpha=0.01)
    
    push!(fs, nonlinear_function.(xs))
end
axes[1].scatter(X_obs, Y_obs)
set_options(axes[1], "x", "y (prob.)", "function samples")

# 予測確率
axes[2].plot(xs, mean(fs), label="prediction")
axes[2].scatter(X_obs, Y_obs, label="data")
set_options(axes[2], "x", "y (prob.)", "prediction", legend=true)

tight_layout()
In [252]:
# 予測対象の点候補
x_list = [-1.0, 0.0, 1.0]

fig_num = length(x_list)
fig, axes = subplots(fig_num, 2, sharey=true, figsize=(12, 4*fig_num))
for (j, x) in enumerate(x_list)
    # パラメータごとに関数を可視化
    for (i, param) in enumerate(param_posterior)
        w1, w2 = param
        nonlinear_function(x) = sig(w1*x + w2)
        axes[j].plot(xs, nonlinear_function.(xs), "g", alpha=0.01)
    end
    
    # 観測データ
    axes[j].scatter(X_obs, Y_obs, label="data")
    
    # 候補点のx座標
    axes[j].plot([x, x], [0, 1], "r--", label="x_p=$(x)")
    
    axes[j].set_xlim(extrema(xs))
    set_options(axes[j], "x", "y (prob.)", "function sample")
    axes[j].legend(loc="lower right")
    
    # 点xにおける関数値(確率値)をヒストグラムとして可視化
    probs = [sig(param[1]*x+param[2]) for param in param_posterior]
    axes[j+fig_num].hist(probs, orientation="horizontal")
    set_options(axes[j+fig_num], "y", "freqency", "function samples";
                            gridy=true)
    axes[j+fig_num].grid(axis="x")
end
In [253]:
#  同時分布p(Y, w|X)
p_joint(X, Y, w1, w2) = prod(pdf.(Bernoulli.(sig.(w1.*X .+ w2)), Y)) * pdf(Normal(μ1, σ1), w1) * pdf(Normal(μ2, σ2), w2)

# μの積分範囲
w_range = range(-30, 30, length=100)

# 数値積分の実行
p_marginal = approx_integration_2D(w_range, p_joint)
p_marginal(X_obs, Y_obs)
Out[253]:
0.2966829565005522
In [254]:
w_posterior = [p_joint(X_obs, Y_obs, w1, w2) for w1 in w_range, w2 in w_range] ./ p_marginal(X_obs, Y_obs)

# 事後分布の描画
fig, axes = subplots(1, 2, figsize=(8,4))
cs = axes[1].contour(w_range, w_range, w_posterior' .+ eps(), cmap="jet")
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "posterior density (contour)")

xgrid = repeat(w_range', length(w_range), 1)
ygrid = repeat(w_range, 1, length(w_range))
axes[2].pcolormesh(xgrid, ygrid, w_posterior', cmap="jet", shading="auto")
set_options(axes[2], "w1", "w2", "posteior density (colored)")

tight_layout()
In [255]:
p_likelihood(xp, yp, w1, w2) = pdf(Bernoulli(sig(w1*xp + w2)), yp)
p_predictive = approx_predictive(w_posterior, w_range, p_likelihood)
Out[255]:
#80 (generic function with 1 method)
In [256]:
# 関数を可視化する範囲
xs = range(-3, 3, length=50)

fig, ax = subplots()
ax.scatter(X_obs, Y_obs, label="data")
ax.plot(xs, p_predictive.(xs, 1), label="probability")
ax.set_xlim([-3, 3])
set_options(ax, "x", "y", "prediction", legend=true)
Out[256]:
PyObject <matplotlib.legend.Legend object at 0x7fe05b723ac0>
In [257]:
using Distributions, PyPlot, LinearAlgebra

# n次元単位行列
# eye(n) = Diagonal{Float64}(I, n)
eye(n) = Diagonal{Float64}(ones(n))

# パラメータ抽出用の関数
unzip(a) = map(a -> getfield.(a,x), fieldnames(eltype(a)))
Out[257]:
unzip (generic function with 1 method)
In [258]:
# 入力データセット
X_obs = [-2, 1, 5]

# 出力データセット
Y_obs = [-2.2, -1.0, 1.5]

# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
set_options(ax, "x", "y", "data (N = $(length(X_obs)))")
Out[258]:
false
In [259]:
# 切片はゼロで固定
w2 = 0

# yに付加されているノイズの標準偏差
σ = 1.0

# 事前分布の平均値と標準偏差
μ1 = 0.0
σ1 = 10.0
Out[259]:
10.0
In [260]:
# 非正規化対数事後分布
ulp(w1) = sum(logpdf.(Normal.(w1*X_obs .+ w2, σ), Y_obs)) + logpdf(Normal(μ1, σ1), w1)

# 分布を表示する範囲
w1s = range(-5, 5, length=100)

fig, axes = subplots(1, 2, figsize=(8,4))

# 非正規化対数事後分布の可視化
axes[1].plot(w1s, ulp.(w1s))
set_options(axes[1], "w1", "log density (unnormalized)", "unnormarized log posterior distribution")

# 非正規化事後分布の可視化
axes[2].plot(w1s, exp.(ulp.(w1s)))
set_options(axes[2], "w1", "density (unnormalized)", "unnormarized posterior distribution")

tight_layout()
In [261]:
# 勾配法(1dim)
function gradient_method_1dim(funct, x_init, η, maxiter)
    ffunct(x) = ForwardDiff.derivative(funct, x)
    x_seq = Array{typeof(x_init), 1}(undef, maxiter)
    x_seq[1] = x_init
    for i in 2:maxiter
        x_seq[i] = x_seq[i-1] + η*ffunct(x_seq[i-1])
    end
    x_seq
end

# 最適化パラメータ
w1_init = 0.0
maxiter = 100
η = 0.01

# 最適化の実施
w1_seq = gradient_method_1dim(ulp, w1_init, η, maxiter)

# 勾配法の過程を可視化
fig, ax = subplots(figsize=(8,4))
ax.plot(w1_seq)
set_options(ax, "iteration", "w1", "w1 sequence")
Out[261]:
false
In [262]:
# 最適化パラメータ
w1_init = 0.0
maxiter = 100
η = 0.01

# 最適化のラッパー関数の定義
function inference_wrapper_gd_1dim(log_joint, params, w1_init, η, maxiter)
    ulp(w1) = log_joint(w1, params...)
    w1_seq = gradient_method_1dim(ulp, w1_init, η, maxiter)
    w1_seq
end

# 対数同時分布
log_joint(w1, X, Y, w2, σ, μ1, σ1) = sum(logpdf.(Normal.(w1*X .+ w2, σ), Y)) + logpdf(Normal(μ1, σ1), w1)
params = (X_obs, Y_obs, w2, σ, μ1, σ1)

# 最適化の実施
w1_seq = inference_wrapper_gd_1dim(log_joint, params, w1_init, η, maxiter)

# 勾配法の過程を可視化
fig, ax = subplots(figsize=(8,4))
ax.plot(w1_seq)
set_options(ax, "iteration", "w1", "w1 sequence")
Out[262]:
false
In [263]:
# 最適化パラメータ
w1_init = 0.0
maxiter = 100
η = 0.01

# グローバル変数を参照する場合

@time gradient_method_1dim(ulp, w1_init, η, maxiter)

@time inference_wrapper_gd_1dim(log_joint, params, w1_init, η, maxiter)
  0.000263 seconds (1.98 k allocations: 76.672 KiB)
  0.000042 seconds (796 allocations: 35.094 KiB)
Out[263]:
100-element Vector{Float64}:
 0.0
 0.109
 0.18528909999999998
 0.23868384109
 0.276054820378891
 0.3022107687831858
 0.3205173170713517
 0.33333007021823907
 0.3422977161457455
 0.34857417153040726
 0.352967062654132
 0.356041647151627
 0.35819354884142374
 ⋮
 0.3632122625791319
 0.3632122625791344
 0.3632122625791362
 0.3632122625791374
 0.3632122625791383
 0.3632122625791389
 0.3632122625791393
 0.36321226257913963
 0.36321226257913986
 0.36321226257913997
 0.3632122625791401
 0.36321226257914013
In [264]:
# 近似分布用の平均を求める
μ_approx = w1_seq[end]

# 分布を表示する範囲
w1s = range(-5, 5, length=100)

fig, axes = subplots(1, 2, figsize=(8,4))

# 非正規化対数事後分布と最適値
axes[1].plot(w1s, ulp.(w1s))
axes[1].plot(μ_approx, ulp(μ_approx), "rx", label="optimal")
set_options(axes[1], "w1", "log density (unnormalized)", "unnormalized log posterior distribution")

# 非正規化事後分布と最適値
axes[2].plot(w1s, exp.(ulp.(w1s)))
axes[2].plot(μ_approx, exp.(ulp(μ_approx)), "rx", label="optimal")
set_options(axes[2], "w1", "density (unnormalized)", "unnormalized posterior distribution")

tight_layout()
In [265]:
# 2階微分から分散を求める
grad(x) = ForwardDiff.derivative(ulp, x)
hessian(x) = ForwardDiff.derivative(grad, x)
σ_approx = sqrt(inv(-hessian(μ_approx)))
Out[265]:
0.1825437644092281
In [266]:
# 可視化する範囲
w1s = range(-0.5, 1.0, length=100)

# 非正規化事後分布の可視化
fig, axes = subplots(1, 2, figsize=(8,4))
axes[1].plot(w1s, exp.(ulp.(w1s)))
set_options(axes[1], "w1", "density (unnormalized)", "unnormalized posterior distributions")

# 得られた近似分布の可視化
axes[2].plot(w1s, pdf.(Normal.(μ_approx, σ_approx), w1s))
set_options(axes[2], "w1", "density", "approximate distribution")

tight_layout()
In [267]:
# 事前分布の設定
σ = 1.0
μ1 = 0.0
μ2 = 0.0
σ1 = 10.0
σ2 = 10.0

# 対数同時分布
log_joint(w, X, Y, σ, μ1, σ1, μ2, σ2) = 
    sum(logpdf.(Normal.(w[1]*X .+ w[2], σ), Y)) + 
    logpdf(Normal(μ1, σ1), w[1]) + 
    logpdf(Normal(μ2, σ2), w[2])

# 非正規化対数事後分布
params = (X_obs, Y_obs, σ, μ1, σ1, μ2, σ2)
ulp(w) = log_joint(w, params...)

# 分布を可視化する範囲
w1s = range(-5, 5, length=100)
w2s = range(-5, 5, length=100)

fig, axes = subplots(1, 2, figsize=(12,4))

# 非正規化対数事後分布の可視化
cs = axes[1].contour(w1s, w2s, [ulp([w1, w2]) for w1 in w1s, w2 in w2s]', cmap="jet")
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "unnormalized log posterior")

# 非正規化事後分布の可視化
cs = axes[2].contour(w1s, w2s, [exp(ulp([w1, w2])) for w1 in w1s, w2 in w2s]', cmap="jet")
axes[2].clabel(cs, inline=true)
set_options(axes[2], "w1", "w2", "unnormalized posterior")

tight_layout()
In [268]:
# 多次元の勾配法
function gradient_method(funct, x_init, η, maxiter)
    x_seq = Array{typeof(x_init[1]), 2}(undef, length(x_init), maxiter)
    g(x) = ForwardDiff.gradient(funct, x)
    x_seq[:, 1] = x_init
    for i in 2:maxiter
        x_seq[:,i] = x_seq[:, i-1] + η*g(x_seq[:, i-1])
    end
    x_seq
end

function inference_wrapper_gd(log_joint, params, w_init, η, maxiter)
    ulp(w) = log_joint(w, params...)
    w_seq = gradient_method(ulp, w_init, η, maxiter)
    w_seq
end

# 最適化パラメータ
w_init = [0.0, 0.0]
maxiter = 1000
η = 0.01

# 最適化の実行
w_seq = inference_wrapper_gd(log_joint, params, w_init, η, maxiter)

# 勾配法の過程を可視化
fig, axes =subplots(2, 1, figsize=(8,8))
axes[1].plot(w_seq[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence")
axes[2].plot(w_seq[2,:])
set_options(axes[2], "iteration", "w2", "w2, sequence")

tight_layout()
In [269]:
# 平均
μ_approx = w_seq[:, end]

# 共分散
hessian(w) = ForwardDiff.hessian(ulp, w)
Σ_approx = inv(-hessian(μ_approx))

fog, axes =subplots(1, 2, figsize=(8,4))

# 非正規化事後分布の可視化
cs = axes[1].contour(w1s, w2s, [exp(ulp([w1, w2])) for w1 in w1s, w2 in w2s]', cmap="jet")
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "unnormalized posterior")

# 近似正規分布の可視化
cx = axes[2].contour(w1s, w2s, [pdf(MvNormal(μ_approx, Σ_approx), [w1, w2]) for w1 in w1s, w2 in w2s]', cmap="jet")
axes[2].clabel(cs, inline=true)
set_options(axes[2], "w1", "w2", "approximate posterior")

tight_layout()
In [270]:
# 長方形の幅
Δ1 = w1s[2] - w1s[1]
Δ2 = w2s[2] - w2s[1]

# 積分近似による予測分布
p_predictive2(x, y) = sum([pdf(Normal(w1*x + w2, σ), y) * 
                                                pdf(MvNormal(μ_approx, Σ_approx), [w1, w2]) * 
                                                Δ1 * Δ2 for w1 in w1s, w2 in w2s])

# 描画範囲
xs = range(-10, 10, length=100)
ys = range(-5, 5, length=100)

# 密度の計算
density_y = p_predictive2.(xs, ys')

fig, axes = subplots(1, 2, sharey=true, figsize=(8,4))

# 等高線図
cs = axes[1].contour(xs, ys, density_y', cmap="jet")
axes[1].clabel(cs, inline=true)
axes[1].scatter(X_obs, Y_obs)
set_options(axes[1], "x", "y", "predictive distribution (contour)")

# カラーメッシュ
xgrid = repeat(xs', length(ys), 1)
ygrid = repeat(ys, 1, length(xs))
axes[2].pcolormesh(xgrid, ygrid, density_y', cmap="jet", shading="auto")
axes[2].plot(X_obs, Y_obs, "ko", label="data")
set_options(axes[2], "x", "y", "predictive distribution (colored)")

tight_layout()
In [271]:
# 入力データ集合
X_obs = [-2, 1, 2]

# 出力データ
Y_obs = Bool.([0, 1, 1])

# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
ax.set_xlim([-3, 3])
set_options(ax, "x", "y", "data (N = $(length(X_obs)))")
Out[271]:
false
In [272]:
# シグモイド関数
sig(x) = 1/(1 + exp(-x))

# 事前分布の設定
μ1 = 0
μ2 = 0
σ1 = 10.0
σ2 = 10.0

# 対数同時分布
log_joint(w, X, Y, μ1, σ1, μ2, σ2) = 
    logpdf(Normal(μ1, σ1), w[1]) + 
    logpdf(Normal(μ2, σ2), w[2]) + 
    sum(logpdf.(Bernoulli.(sig.(w[1]*X .+ w[2])), Y))
params = (X_obs, Y_obs, μ1, σ1, μ2, σ2)

# 非正規化対数事後分布
ulp(w) = log_joint(w, params...)
Out[272]:
ulp (generic function with 1 method)
In [273]:
# 最適化パラメータ
w_init = [0.0, 0.0]
maxiter = 2000
η = 0.1

# 最適化の実行
w_seq = inference_wrapper_gd(log_joint, params, w_init, η, maxiter)

# 勾配法の過程を可視化
fig, axes = subplots(2, 1, figsize=(8,8))
axes[1].plot(w_seq[1,:])
set_options(axes[1], "iteration", "w1", "w1 sequence")
axes[2].plot(w_seq[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence")

tight_layout()
In [274]:
# 平均
μ_approx = w_seq[:, end]

# 共分散
hessian(x) = ForwardDiff.hessian(ulp, x)
Σ_approx = inv(-hessian(μ_approx))
Out[274]:
2×2 Matrix{Float64}:
 18.6137   -2.69261
 -2.69261  30.6517
In [275]:
w1s = range(-10, 30, length=100)
w2s = range(-20, 20, length=100)

fig, axes = subplots(1, 2, figsize=(8,4))
cs = axes[1].contour(w1s, w2s, [exp(ulp([w1, w2])) + eps() for w1 in w1s, w2 in w2s]', cmap="jet")
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "unnormalized posterior")

cs = axes[2].contour(w1s, w2s, [pdf(MvNormal(μ_approx, Σ_approx), [w1, w2]) for w1 in w1s, w2 in w2s]', cmap="jet")
axes[2].clabel(cs, inline=true)
set_options(axes[2], "w1", "w2", "approximate posterior")

tight_layout()
In [276]:
# 近似分布から候補のサンプルを100個抽出
W = rand(MvNormal(μ_approx, Σ_approx), 100)

fig, ax = subplots()
for i in 1:size(W, 2)
    w1, w2 = W[:, i]
    funct(x) = sig(w1*x + w2)
    ax.plot(xs, funct.(xs), "g", alpha=0.1)
end
ax.scatter(X_obs, Y_obs)
ax.set_xlim(extrema(xs))
set_options(ax, "x", "y", "prediction sample from approximate posterior")
Out[276]:
false
In [277]:
# 長方形の幅
Δ1 = w1s[2] - w1s[1]
Δ2 = w2s[2] - w2s[1]

# 積分近似による予測分布
p_predictive2(x, y) = sum([pdf(Bernoulli(sig(w1*x + w2)), y) * 
        pdf(MvNormal(μ_approx, Σ_approx), [w1, w2]) * 
        Δ1 * Δ2 for w1 in w1s, w2 in w2s])

# 関数を可視化する範囲
xs = range(-3, 3, length=50)

fig, ax = subplots()
ax.scatter(X_obs, Y_obs, label="data")
ax.plot(xs, p_predictive2.(xs, 1), label="probability")
ax.set_xlim([-3, 3])
set_options(ax, "x", "y", "predection", legend=true)
Out[277]:
PyObject <matplotlib.legend.Legend object at 0x7fe03ef0ae20>
In [278]:
# ガウス提案分布によるメトロポリスヘイスティング
function GaussianMH(log_p_tilde, μ0; maxiter::Int=100_000, σ::Float64=1.0)
    # サンプルを格納する配列
    D = length(μ0)
    μ_samples = Array{typeof(μ0[1]), 2}(undef, D, maxiter)
    
    #  初期サンプル
    μ_samples[:, 1] = μ0
    
    # 受容されたサンプルの数
    num_accepted = 1
    
    for i in 2:maxiter
        # 提案分布q(多変量正規分布)に従い、次のサンプルの候補を抽出
        μ_tmp = rand(MvNormal(μ_samples[:, i-1], σ*eye(D)))
        
        # 比率r(の対数)を計算
        log_r = (log_p_tilde(μ_tmp) + logpdf(MvNormal(μ_tmp, σ*eye(D)), μ_samples[:, i-1])) - 
                        (log_p_tilde(μ_samples[:, i-1]) + logpdf(MvNormal(μ_samples[:, i-1], σ*eye(D)), μ_tmp))
        
        # 確率rでサンプルを受容する
        is_accepted = min(1, exp(log_r)) > rand()
        new_sample = is_accepted ? μ_tmp : μ_samples[:, i-1]
        
        # 新しいサンプルを格納
        μ_samples[:, i] = new_sample
        
        # 受容された場合、合計をプラスする
        num_accepted += is_accepted
    end
    
    μ_samples, num_accepted
end
Out[278]:
GaussianMH (generic function with 1 method)
In [279]:
# ハミルトニアンモンテカルロ法
function HMC(log_p_tilde, μ0; maxiter::Int=100_000, L::Int=100, ϵ::Float64=1e-1)
    # leapflogによる値の更新
    function leapflog(grad, p_in, μ_in, L, ϵ)
        μ = μ_in
        p = p_in + 0.5*ϵ*grad(μ)
        for l in 1:L-1
            μ += ϵ*p
            p += ϵ*grad(μ)
        end
        μ += ϵ*p
        p += 0.5*ϵ*grad(μ)
        p, μ
    end
    
    # 非正規化対数事後分布の勾配関数の計算
    grad(μ) = ForwardDiff.gradient(log_p_tilde, μ)
    
    # サンプルを格納する配列
    D = length(μ0)
    μ_samples = Array{typeof(μ0[1]), 2}(undef, D, maxiter)
    
    # 初期サンプル
    μ_samples[:, 1] = μ0
    
    # 受容されたサンプルの数
    num_accepted = 1
    
    for i in 2:maxiter
        # 運動量pの生成
        p_in = randn(size(μ0))
        
        # リープフロッグ法
        p_out, μ_out = leapflog(grad, p_in, μ_samples[:, i-1], L, ϵ)
        
        # 比率r(の対数)を計算
        μ_in = μ_samples[:, i-1]
        log_r = (log_p_tilde(μ_out) + logpdf(MvNormal(zeros(D), eye(D)), vec(p_out))) -
                        (log_p_tilde(μ_in) + logpdf(MvNormal(zeros(D), eye(D)), vec(p_in)))
        
        # 確率rでサンプルを受容する
        is_accepted = min(1, exp(log_r)) > rand()
        new_sample = is_accepted ? μ_out : μ_in
        
        # 新しいサンプルを受容する
        μ_samples[:, i] = new_sample
        
        # 受容された場合、合計をプラスする
        num_accepted += is_accepted
    end
    
    μ_samples, num_accepted
end
Out[279]:
HMC (generic function with 1 method)
In [280]:
function inference_wrapper_GMH(log_joint, params, w_init; maxiter::Int=100_000, σ::Float64=1.0)
    ulp(w) = log_joint(w, params...)
    GaussianMH(ulp, w_init; maxiter=maxiter, σ=σ)
end

function inference_wrapper_HMC(log_joint, params, w_init; maxiter::Int=100_000, L::Int=100, ϵ::Float64=1e-1)
    ulp(w) = log_joint(w, params...)
    HMC(ulp, w_init; maxiter=maxiter, L=L, ϵ=ϵ)
end
Out[280]:
inference_wrapper_HMC (generic function with 1 method)
In [281]:
# 入力データセット
X_obs = [-2, 1, 5]

# 出力データセット
Y_obs = [-2.2, -1.0, 1.5]

# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
set_options(ax, "x", "y", "data (N = $(length(X_obs)))")
Out[281]:
false
In [282]:
# 対数同時分布
log_joint(w, X, Y, σ, μ1, σ1, μ2, σ2) = 
    sum(logpdf.(Normal.(w[1]*X .+ w[2], σ), Y)) + 
    logpdf(Normal(μ1, σ1), w[1]) + 
    logpdf(Normal(μ2, σ2), w[2])
params = (X_obs, Y_obs, σ, μ1, σ1, μ2, σ2)

# 非正規化対数事後分布
ulp(w) = log_joint(w, params...)
Out[282]:
ulp (generic function with 1 method)
In [283]:
# 初期値
w_init = randn(2)

# サンプリング
maxiter = 300
param_posterior_GMH, num_accepted_GMH = 
    inference_wrapper_GMH(log_joint, params, w_init, maxiter=maxiter, σ=1.0)
param_posterior_HMC, num_accepted_HMC = 
    inference_wrapper_HMC(log_joint, params, w_init, maxiter=maxiter, L=10, ϵ=1e-1)

# サンプリングの過程を可視化(Gaussian MH)
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior_GMH[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence (GMH)")
axes[2].plot(param_posterior_GMH[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence (GMH)")
tight_layout()
println("acceptance rate (GMH) = $(num_accepted_GMH/maxiter)")

# サンプリングの過程を可視化(HMC)
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior_HMC[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence (HMC)")
axes[2].plot(param_posterior_HMC[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence (HMC)")
tight_layout()
println("acceptance rate (HMC) = $(num_accepted_HMC/maxiter)")
acceptance rate (GMH) = 0.16
acceptance rate (HMC) = 0.9933333333333333
In [284]:
# 事後分布を可視化する範囲
w1s = range(-5, 5, length=100)
w2s = range(-5, 5, length=100)

fig, axes =subplots(1, 3, sharex=true, sharey=true, figsize=(12,4))
cs = axes[1].contour(w1s, w2s, [exp(ulp([w1, w2])) + eps() for w1 in w1s, w2 in w2s]', cmap="jet")

# 非正規化事後分布
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "unnormalized posterior")

# Gaussian MHで得られた事後分布からのサンプル
axes[2].scatter(param_posterior_GMH[1, :], param_posterior_GMH[2, :], alpha=10/maxiter)
set_options(axes[2], "w1", "w2", "samples from posterior (GMH)")

# HMCで得られた事後分布からのサンプル
axes[3].scatter(param_posterior_HMC[1, :], param_posterior_HMC[2, :], alpha=10/maxiter)
set_options(axes[3], "w1", "w2", "samples from posterior (HMC)")
Out[284]:
false
In [285]:
xs = range(-10, 10, length=100)

fig, axes = subplots(1, 2, figsize=(8,4))

# Gaussian MH
for i in 1:size(param_posterior_GMH, 2)
    w1, w2 = param_posterior_GMH[:, i]
    funct(x) = w1*x + w2
    axes[1].plot(xs, funct.(xs), "g", alpha=10/maxiter)
end
axes[1].scatter(X_obs, Y_obs)
axes[1].set_xlim(extrema(xs))
set_options(axes[1], "x", "y", "predictive distribution (GMH)")

# HMC
for i in 1:size(param_posterior_HMC, 2)
    w1, w2 = param_posterior_HMC[:, i]
    funct(x) = w1*x + w2
    axes[2].plot(xs, funct.(xs), "g", alpha=10/maxiter)
end
axes[2].scatter(X_obs, Y_obs)
axes[2].set_xlim(extrema(xs))
set_options(axes[2], "x", "y", "predictive distribution (HMC)")
Out[285]:
false
In [286]:
# 入力データセット
X_obs = [-2, 1, 2]

# 出力データセット
Y_obs = Bool.([0, 1, 1])

# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
set_options(ax, "x", "y", "data = (N = $(length(X_obs)))")
Out[286]:
false
In [287]:
# シグモイド関数
sig(x) = 1/(1 + exp(-x))

# 事前分布の設定
μ1 = 0
μ2 = 0
σ1 = 10.0
σ2 = 10.0

# 対数同時分布
log_joint(w, X, Y, μ1, σ1, μ2, σ2) = 
    logpdf(Normal(μ1, σ1), w[1]) + 
    logpdf(Normal(μ2, σ2), w[2]) + 
    sum(logpdf.(Bernoulli.(sig.(w[1]*X .+ w[2])), Y))
params = (X_obs, Y_obs, μ1, σ1, μ2, σ2)

# 非正規化対数事後分布
ulp(w) = log_joint(w, params...)
Out[287]:
ulp (generic function with 1 method)
In [288]:
# 初期値
w_init = randn(2)

# サンプリング
maxiter = 300
param_posterior_GMH, num_accepted_GMH = 
    inference_wrapper_GMH(log_joint, params, w_init, maxiter=maxiter, σ=1.0)
param_posterior_HMC, num_accepted_HMC = 
    inference_wrapper_HMC(log_joint, params, w_init, maxiter=maxiter, L=10, ϵ=1e-1)

# サンプリングの過程を可視化(Gaussian MH)
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior_GMH[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence (GMH)")
axes[2].plot(param_posterior_GMH[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence (GMH)")
tight_layout()
println("acceptance rate (GMH) = $(num_accepted_GMH/maxiter)")

# サンプリングの過程を可視化(HMC)
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior_HMC[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence (HMC)")
axes[2].plot(param_posterior_HMC[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence (HMC)")
tight_layout()
println("acceptance rate (HMC) = $(num_accepted_HMC/maxiter)")
acceptance rate (GMH) = 0.9033333333333333
acceptance rate (HMC) = 1.0
In [289]:
# 事後分布を可視化する範囲
w1s = range(-10, 30, length=100)
w2s = range(-20, 20, length=100)

fig, axes =subplots(1, 3, sharex=true, sharey=true, figsize=(12,4))
cs = axes[1].contour(w1s, w2s, [exp(ulp([w1, w2])) + eps() for w1 in w1s, w2 in w2s]', cmap="jet")

# 非正規化事後分布
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "unnormalized posterior")

# Gaussian MHで得られた事後分布からのサンプル
axes[2].scatter(param_posterior_GMH[1, :], param_posterior_GMH[2, :], alpha=10/maxiter)
set_options(axes[2], "w1", "w2", "samples from posterior (GMH)")

# HMCで得られた事後分布からのサンプル
axes[3].scatter(param_posterior_HMC[1, :], param_posterior_HMC[2, :], alpha=10/maxiter)
set_options(axes[3], "w1", "w2", "samples from posterior (HMC)")
Out[289]:
false
In [290]:
param_posterior_GMH, num_accepted_GMH = 
    GaussianMH(ulp, w_init, maxiter=maxiter, σ=10.0)
param_posterior_HMC, num_accepted_HMC = 
    HMC(ulp, w_init, maxiter=maxiter, L=10, ϵ=1e-0)

# サンプリングの過程を可視化(Gaussian MH)
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior_GMH[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence (GMH)")
axes[2].plot(param_posterior_GMH[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence (GMH)")
tight_layout()
println("acceptance rate (GMH) = $(num_accepted_GMH/maxiter)")

# サンプリングの過程を可視化(HMC)
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior_HMC[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence (HMC)")
axes[2].plot(param_posterior_HMC[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence (HMC)")
tight_layout()
println("acceptance rate (HMC) = $(num_accepted_HMC/maxiter)")
acceptance rate (GMH) = 0.6466666666666666
acceptance rate (HMC) = 0.9766666666666667
In [291]:
# 事後分布を可視化する範囲
w1s = range(-10, 30, length=100)
w2s = range(-20, 20, length=100)

fig, axes =subplots(1, 3, sharex=true, sharey=true, figsize=(12,4))
cs = axes[1].contour(w1s, w2s, [exp(ulp([w1, w2])) + eps() for w1 in w1s, w2 in w2s]', cmap="jet")

# 非正規化事後分布
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "unnormalized posterior")

# Gaussian MHで得られた事後分布からのサンプル
axes[2].scatter(param_posterior_GMH[1, :], param_posterior_GMH[2, :], alpha=10/maxiter)
set_options(axes[2], "w1", "w2", "samples from posterior (GMH)")

# HMCで得られた事後分布からのサンプル
axes[3].scatter(param_posterior_HMC[1, :], param_posterior_HMC[2, :], alpha=10/maxiter)
set_options(axes[3], "w1", "w2", "samples from posterior (HMC)")
Out[291]:
false
In [292]:
# 関数を可視化する範囲
xs = range(-3, 3, length=100)

fig, axes = subplots(2, 2, figsize=(12,8))

# Gaussian MHによるサンプル
fs =[]
for i in 1:size(param_posterior_GMH, 2)
    w1, w2 = param_posterior_GMH[:, i]
    funct(x) = sig(w1*x + w2)
    push!(fs, funct.(xs))
    axes[1,1].plot(xs, funct.(xs), "g", alpha=10/maxiter)
end
axes[1,1].scatter(X_obs, Y_obs)
axes[1,1].set_xlim(extrema(xs))
set_options(axes[1,1], "x", "y", "predictive distribution (GMH)")

# GMHによる予測の平均
axes[1,2].plot(xs, mean(fs), label="prediction")
axes[1,2].scatter(X_obs, Y_obs, label="data")
set_options(axes[1,2], "x", "y (prob.)", "prediction (GMH)", legend=true)

# HMCによるサンプル
fs =[]
for i in 1:size(param_posterior_HMC, 2)
    w1, w2 = param_posterior_HMC[:, i]
    funct(x) = sig(w1*x + w2)
    push!(fs, funct.(xs))
    axes[2,1].plot(xs, funct.(xs), "g", alpha=10/maxiter)
end
axes[2,1].scatter(X_obs, Y_obs)
axes[2,1].set_xlim(extrema(xs))
set_options(axes[2,1], "x", "y", "predictive distribution (HMC)")

# HMCによる予測の平均
axes[2,2].plot(xs, mean(fs), label="prediction")
axes[2,2].scatter(X_obs, Y_obs, label="data")
set_options(axes[2,2], "x", "y (prob.)", "prediction (HMC)", legend=true)
Out[292]:
PyObject <matplotlib.legend.Legend object at 0x7fe055642130>
In [293]:
# 入力データセット
X_obs = [-2, -1.5, 0.5, 0.7, 1.2]

# 出力データセット
Y_obs = [0, 0, 2, 1, 8]

# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
ax.set_xlim([-2, 2])
set_options(ax, "x", "y", "data = (N = $(length(X_obs)))")
Out[293]:
false
In [294]:
# 対数同時分布
log_joint(w, X, Y) = 
    sum(logpdf.(Poisson.(exp.(w[1]*X .+ w[2])), Y)) +
    logpdf(Normal(μ1, σ1), w1) + 
    logpdf(Normal(μ2, σ2), w2)
params = (X_obs, Y_obs)

# 非正規化対数事後分布
ulp(w) = log_joint(w, params...)

# 初期値
w_init = randn(2)

# サンプルサイズ
maxiter = 300

# ハミルトニアンモンテカルロ法によるサンプリング
param_posterior, num_accepted = inference_wrapper_HMC(log_joint, params, w_init, maxiter=maxiter)

# トレースプロット
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence")
axes[2].plot(param_posterior[2, :])
set_options(axes[2], "iteration", "w2", "w2 sequence")
tight_layout()
println("acceptance rate = $(num_accepted/maxiter)")

# HMCで得られた事後分布からのサンプル
fig, ax = subplots()
ax.scatter(param_posterior[1, :], param_posterior[2, :], alpha=0.1)
set_options(ax, "w1", "w2", "sample from posterior")
tight_layout()
acceptance rate = 0.97
In [295]:
# 関数を可視化する範囲
xs = range(-2, 2, length=100)

fig, ax = subplots()
for i in 1:size(param_posterior, 2)
    w1, w2 = param_posterior[:, i]
    
    # 指数関数を可視化
    funct(x) = exp(w1*x + w2)
    ax.plot(xs, funct.(xs), "g", alpha=0.1)
end
ax.plot(X_obs, Y_obs, "ko")
ax.set_xlim(extrema(xs))
ax.set_ylim([-1, 15])
set_options(ax, "x", "y", "predictive distribution")

tight_layout()
In [296]:
# 学習データ
X_obs = []
Y_obs = []

push!(X_obs, [0.3, 0.4])
push!(Y_obs, [4.0, 3.7])

push!(X_obs, [0.2, 0.4, 0.9])
push!(Y_obs, [6.0, 7.2, 9.4])

push!(X_obs, [0.6, 0.8, 0.9])
push!(Y_obs, [6.0, 6.9, 7.8])

# 散布図の可視化
fig, ax = subplots()
ax.plot(X_obs[1], Y_obs[1], "or", label="class1")
ax.plot(X_obs[2], Y_obs[2], "xg", label="class2")
ax.plot(X_obs[3], Y_obs[3], "^b", label="class3")
set_options(ax, "x", "y", "data", legend=true)
Out[296]:
PyObject <matplotlib.legend.Legend object at 0x7fe03ea39e80>
In [297]:
function linear_fit(Y, X)
    N = length(Y)
    w1 = sum((Y .- mean(Y)) .* X) / sum((X .- mean(X)) .* X)
    w2 = mean(Y) - w1*mean(X)
    w1, w2
end

# まとめて回帰
w1, w2 = linear_fit(vcat(Y_obs...), vcat(X_obs...))

# 個別に回帰
w1s = []
w2s = []
for i in 1:3
    w1_tmp, w2_tmp = linear_fit(Y_obs[i], X_obs[i])
    push!(w1s, w1_tmp)
    push!(w2s, w2_tmp)
end

# 関数を可視化する範囲
xs = range(0, 1, length=100)

fig, axes = subplots(1, 2, figsize=(8,4))

# 単一の回帰
axes[1].plot(xs, w1.*xs .+ w2, "-k")

# 個別の回帰
axes[2].plot(xs, w1s[1].*xs .+ w2s[1], "-r")
axes[2].plot(xs, w1s[2].*xs .+ w2s[2], "-g")
axes[2].plot(xs, w1s[3].*xs .+ w2s[3], "-b")

# データの可視化
for ax in axes
    ax.plot(X_obs[1], Y_obs[1], "or", label="class1")
    ax.plot(X_obs[2], Y_obs[2], "xg", label="class2")
    ax.plot(X_obs[3], Y_obs[3], "^b", label="class3")
end

set_options(axes[1], "x", "y", "(a) single regression", legend=true)
set_options(axes[2], "x", "y", "(b) multiple regressions", legend=true)

tight_layout()
In [298]:
# 対数同時分布の設計
@views hyper_prior(w) = logpdf(Normal(0, 10.0), w[1]) + logpdf(Normal(0, 10.0), w[2])
@views prior(w) = sum(logpdf.(Normal.(w[1], 1.0), w[3:5])) + sum(logpdf.(Normal.(w[2], 1.0), w[6:8]))
@views log_likelihood(Y, X, w) = sum([sum(logpdf.(Normal.(Y[i], 1.0), w[2+i].*X[i] .+ w[2+i+3])) for i in 1:3])
log_joint(w, X, Y) = hyper_prior(w) + prior(w) + log_likelihood(Y, X, w)
params = (X_obs, Y_obs)
ulp(w) = log_joint(w, params...)
Out[298]:
ulp (generic function with 1 method)
In [299]:
# HMCによるサンプリング
maxiter = 1000
w_init = randn(8)
param_posterior, num_accepted = inference_wrapper_HMC(log_joint, params, w_init, maxiter=maxiter, L=100, ϵ=1e-2)

# 予測分布の可視化
fig, axes = subplots(1, 3, sharey=true, figsize=(12, 4))

for i in 1:3
    for j in 1:size(param_posterior, 2)
        w1, w2 = param_posterior[[2+i, 2+i+3], j]
        axes[i].plot(xs, w1.*xs .+ w2, "c-", alpha=0.01)
    end
    set_options(axes[i], "x", "y", "class $(i)")
end
axes[1].plot(X_obs[1], Y_obs[1], "or")
axes[2].plot(X_obs[2], Y_obs[2], "xg")
axes[3].plot(X_obs[3], Y_obs[3], "^b")

tight_layout()
In [300]:
X_obs = []
Y_obs = []

push!(X_obs, [0.1, 0.3, 0.4, 0.5, 0.6, 0.9])
push!(Y_obs, [4.0, 4.0, 3.7, 3.8, 3.9, 3.7])

push!(X_obs, [0.2, 0.4, 0.9])
push!(Y_obs, [6.0, 7.2, 9.4])

push!(X_obs, [0.6, 0.8, 0.9])
push!(Y_obs, [6.0, 6.9, 7.8])

# 対数同時分布の設計
@views hyper_prior(w) = logpdf(Normal(0, 10.0), w[1]) + logpdf(Normal(0, 10.0), w[2])
@views prior(w) = sum(logpdf.(Normal.(w[1], 1.0), w[3:5])) + sum(logpdf.(Normal.(w[2], 1.0), w[6:8]))
@views log_likelihood(Y, X, w) = sum([sum(logpdf.(Normal.(Y[i], 1.0), w[2+i].*X[i] .+ w[2+i+3])) for i in 1:3])
log_joint(w, X, Y) = hyper_prior(w) + prior(w) + log_likelihood(Y, X, w)
params = (X_obs, Y_obs)
ulp(w) = log_joint(w, params...)

# HMCによるサンプリング
maxiter = 1000
w_init = randn(8)
param_posterior, num_accepted = inference_wrapper_HMC(log_joint, params, w_init, maxiter=maxiter, L=100, ϵ=1e-2)

# 予測分布の可視化
fig, axes = subplots(1, 3, sharey=true, figsize=(12, 4))

for i in 1:3
    for j in 1:size(param_posterior, 2)
        w1, w2 = param_posterior[[2+i, 2+i+3], j]
        axes[i].plot(xs, w1.*xs .+ w2, "c-", alpha=0.01)
    end
    set_options(axes[i], "x", "y", "class $(i)")
end
axes[1].plot(X_obs[1], Y_obs[1], "or")
axes[2].plot(X_obs[2], Y_obs[2], "xg")
axes[3].plot(X_obs[3], Y_obs[3], "^b")

tight_layout()
In [301]:
# 系列の長さ
N = 20

# 観測データの次元
D = 2

# 系列データ(#=, =#は複数業をつなげるために挿入)
Y_obs = 
    [1.9 0.2 0.1 1.4 0.3 1.3 1.6 1.5 1.6 2.4 #=
=# 1.7 3.6 2.8 1.6 3.0 2.8 5.1 5.2 6.0 6.4;
      0.1 0.2 0.9 1.5 4.0 5.0 6.3 5.8 6.4 7.5 #=
=# 6.7 7.6 8.7 8.2 8.5 9.6 8.4 8.4 8.4 9.0]

# 2次元に系列データを可視化
fig, ax = subplots(figsize=(6,6))
ax.plot(Y_obs[1, :], Y_obs[2, :], "kx--")
ax.text(Y_obs[1, 1], Y_obs[2, 1], "start", color="r")
ax.text(Y_obs[1, end], Y_obs[2, end], "goal", color="r")
set_options(ax, "y1", "y2", "2dim data")
Out[301]:
false
In [302]:
# 初期状態に与えるノイズ量
σ1 = 100.0

# 状態の遷移に仮定するノイズ量
σ_x = 1.0

# 観測に仮定するノイズ量
σ_y = 1.0

# 状態の遷移系列に関する対数密度
@views transition(X, σ1, σ_x, D, N) = 
    logpdf(MvNormal(zeros(D), σ1 * eye(D)), X[:, 1]) + 
    sum([logpdf(MvNormal(X[:, n-1], σ_x * eye(D)), X[:, n]) for n in 2:N])

# 観測データに関する対数密度
@views observation(X, Y, σ_y, D, N) = 
    sum([logpdf(MvNormal(X[:, n], σ_y * eye(D)), Y[:, n]) for n in 1:N])

# 対数同時分布
log_joint_tmp(X, Y, σ1, σ_x, σ_y, D, N) = 
    transition(X, σ1, σ_x, D, N) + 
    observation(X, Y, σ_y, D, N)

# DN次元ベクトルを入力とする関数にする
log_joint(X_vec, Y, σ1, σ_x, σ_y, D, N) = 
    log_joint_tmp(reshape(X_vec, D, N), Y, σ1, σ_x, σ_y, D, N)
params = (Y_obs, σ1, σ_x, σ_y, D, N)

# 非正規化対数事後分布
ulp(X_vec) = log_joint(X_vec, params...)
Out[302]:
ulp (generic function with 1 method)
In [303]:
#  初期値
X_init = randn(D * N)

# サンプルサイズ
maxiter = 1000

# HMCの実行
samples, num_accepted = inference_wrapper_HMC(log_joint, params, X_init, maxiter=maxiter)

println("acceptance rate = $(num_accepted/maxiter)")
acceptance rate = 0.995
In [304]:
fig, ax = subplots(figsize=(6,6))
for i in 1:maxiter
    # i番目のサンプルの状態変数
    X = reshape(samples[:, i], D, N)
    ax.plot(X[1, :], X[2, :], "go--", alpha=10.0/maxiter)
end

# 観測データの可視化
ax.plot(Y_obs[1, :], Y_obs[2, :], "kx--", label="observation (Y)")

# 状態変数の平均
mean_trace = reshape(mean(samples, dims=2), D, N)
ax.plot(mean_trace[1, :], mean_trace[2, :], "ro--", label="estimated position (X)")

set_options(ax, "y1", "y2", "2dim data", legend=true)
Out[304]:
PyObject <matplotlib.legend.Legend object at 0x7fe041bbb8b0>
In [305]:
# 観測データ数
N = 20

# 入力値
Z_obs = [10, 10, 10, 10, 10, 10, 10, 10, 10, 15, 15, 15, 15, 15, 15, 15, 8, 8, 8, 8]

# 出力値
Y_obs = [67, 64, 60, 60, 57, 54, 51, 51, 49, 63, 62, 62, 58, 57, 53, 51, 24, 22, 23, 19]

# データの可視化
fig, ax = subplots()
ax.scatter(Z_obs, Y_obs)
set_options(ax, "z", "y", "data (scatter)")
Out[305]:
false
In [306]:
fig, axes = subplots(2, 1, figsize=(8,6))

# 出力値
axes[1].plot(Y_obs)
set_options(axes[1], "time", "y", "time series (Y_obs)")

# 入力値
axes[2].plot(Z_obs)
set_options(axes[2], "time", "z", "time series (Z_obs)")
Out[306]:
false
In [307]:
#  初期状態に与えるノイズ量
σ1 = 10.0

# 状態の遷移に仮定するノイズ量
σ_x = 1.0

# 観測に仮定するノイズ量
σ_y = 0.5

# パラメータに仮定するノイズ量
σ_w = 100.0

# パラメータの事前分布
prior(w, σ_w) = logpdf(MvNormal(zeros(2), σ_w*eye(2)), w)

# 状態の遷移系列に関する対数密度
@views transition(X, σ1, σ_x, N) = 
    logpdf(Normal(0, σ1), X[1]) + 
    sum(logpdf.(Normal.(X[1: N-1], σ_x), X[2:N]))

# 観測データに関する対数密度
@views observation(X, Y, Z, w, σ_y) = 
    sum(logpdf.(Normal.(w[1].*Z .+ w[2] + X, σ_y), Y))

# 対数同時分布
log_joint_tmp(X, w, Y, Z, σ1, σ_x, σ_y, σ_w, N) = 
    transition(X, σ1, σ_x, N) + 
    observation(X, Y, Z, w, σ_y) + 
    prior(w, σ_w)
@views log_joint(X_vec, Y, Z, σ1, σ_x, σ_y, σ_w, N) = 
    log_joint_tmp(X_vec[1:N], X_vec[N+1:N+2], Y, Z, σ1, σ_x, σ_y, σ_w, N)
# log_joint(X_vec, Y, Z, σ1, σ_x, σ_y, σ_w, N) = 
#     transition(X_vec[1:N], σ1, σ_x, N) + 
#     observation(X_vec[1:N], Y, Z, X_vec[N+1:N+2]) + 
#     prior(X_vec[N+1,N+2], σ_w)
params = (Y_obs, Z_obs, σ1, σ_x, σ_y, σ_w, N)

# HMCの実行
x_init = randn(N+2)
maxiter = 1000
samples, num_accepted = inference_wrapper_HMC(log_joint, params, x_init, maxiter=maxiter, L=100, ϵ=1e-2)
println("acceptance rate $(num_accepted/maxiter)")
acceptance rate 0.964
In [308]:
# 推定結果の可視化
fig, axes = subplots(5, 1, figsize=(8,15))

# 出力値
axes[1].plot(Y_obs)
set_options(axes[1], "time", "y", "output data (Y)")

# 入力値
axes[2].plot(Z_obs)
set_options(axes[2], "time", "z", "input data (Z)")

# 回帰によって説明される成分
for i in 1:maxiter
    w1, w2 = samples[N+1:N+2, i]
    axes[3].plot(w1*Z_obs .+ w2, "g-", alpha=10/maxiter)
end
set_options(axes[3], "time", "w1x + w2", "regression")

# 状態変数によって説明される成分
for i in 1:maxiter
    X = samples[1:N,i]
    axes[4].plot(X, "g-", alpha=10/maxiter)
end
set_options(axes[4], "time", "x", "time series")

# 回帰と状態変数の和
for i in 1:maxiter
    w1, w2 = samples[N+1:N+2, i]
    X = samples[1:N, i]
    axes[5].plot(w1*Z_obs .+ w2 + X, "g-", alpha=10/maxiter)
end
axes[5].plot(Y_obs, "k", label="Y_obs")
set_options(axes[5], "time", "w1x + w2 + x", "regression + time series", legend=true)

tight_layout()